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ABSTRACT 

According  to  the  linearized  solutions  for  thermal 

blooming,  the  density  perturbations  become  infinite  (i.e. 

"catastrophic"  defocusing)  as  the  Mach  number  approaches 

unity.  However,  the  nonlinearities  in  the  transonic  equations 

cutoff  the  trend  to  infinity,  and  the  values  of  the  flow 

perturbation  quantities  are  finite.  The  nonlinear  equations 

with  heat  addition  are  transformed  into  simple  linear 

algebraic  equations  through  the  specification  of  the 

streamline  geometry  in  the  heat  release  region.  At  a  Mach 

number  of  unity,  streamtube  area  variation  was  found  to  be 

directly  proportional  to  the  change  in  total  temperature. 

A  steady,  two-dimensional  mixed  flow  solution  has  been  found 

for  the  transonic  thermal  blooming  problem.  The  solution 

for  the  density  perturbations  within  a  laser  beam  at  a  Mach 

number  of  precisely  unity  is  given.  For  a  Gaussian  beam 
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with  an  intensity  of  3.333x10  Watts/ra  and  an  atmospheric 
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absorption  of  8 . 0x10  cm  the  maximum  fractional  density 
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perturbation  is  1.028x10  .  The  transonic  thermal  blooming 

problem  does  not  pose  as  serious  a  problem  as  previously 
anticipated. 
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English 

A  area,  square  centimeters 

a  speed  of  sound,  meters/second 

b  function,  dimensionless 

C  line  contour  around  integration  region  R 

c  speed  of  light,  meters/second 

Cp  specific  heat  at  constant  pressure, 

Joules/kilogram  mass  -  degree  Kelvin 

Cv  specific  heat  at  constant  volume. 

Joules /kilogram  mass  -  degree  Kelvin 

f  function,  dimensionless 

G  Green’s  function,  dimensionless;  also  specified 

function,  dimensionless 

H  Heat  quantity,  dimensionless 

h  specific  enthalpy.  Joules/kilogram  mass 

h(x,y)  heat  addition  function,  Watts/cubic  centimeters 

I  laser  beam  intensity,  Watts/square  centimeter? 

also  integral  function  representation 

I  (x)  unit  function?  zero  for  x  <  0  and  unity  for  x  >_  0 

k  Glads tone-Dale  constant,  cubic  meters/kilogram  mass 

L  scaling  factor,  dimensionless;  length,  meters? 

and  functional  representation  of  normalized  velocity, 
dimensionless 

L  characteristic  length  (subscript  1  for  x-direction 

and  2  for  y-direction) ,  centimeters 

M  Mach  number,  dimensionless 

N  Number  of  waves,  dimensionless 

n  normal  coordinate,  centimeters;  also  index  of 

refraction,  dimensionless 
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p  static  pressure,  Newtons/square  meter 

Pr  Prandtl  number,  dimensionless 

Q  heat  source,  Watts/cubic  centimeter 

q  heat  source,  Watts /square  centimeter?  also, 

strength  of  line  heat  source,  Watts/centimeter 

R  universal  gas  constant,  Joules/kilogram  mass  - 

degrees  Kelvin  ;  radius  of  curvature,  centimeters 
and  region  of  integration 

r  distance,  centimeters 

Re  Reynolds  number,  dimensionless 

S  line  contour  around  shock  waves?  also,  distance 

along  a  characteristic,  centimeters 

s  specific  entropy.  Joules/kilogram  mass  -  degree 

Kelvin  ;  also  streamline  coordinate 

T  static  temperature,  degrees  Kelvin 

t  maximum  streamtube  thickness,  centimeters 

U  velocity  vector,  meters/second 

U,u  speed  (x-direction) ,  meters/second 

v  speed  (y-direction) ,  meters/second 

w  laser  beam  spot  size,  centimeters? 

also  speed,  meters/second 

x  coordinate,  centimeters 

y  coordinate,  centimeters 

Z  streamtube  shape,  centimeters 


t 


Greek 

a  atmospheric  absorption,  1/centimeter s 

2  1/2 

3  variable,  dimensionless?  (1-M  )  '  for  subsonic 

flow,  dimensionless? 

2  1/2 

and  (M  -  1)  for  supersonic  flow,  dimensionless 
y  ratio  of  specific  heats,  dimensionless 
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A 

5 

e 

n 

e 

K 

X 

u 

V 

5 

p 

E 

a 

T 

<5 

4> 

t 

Q 


small  change,  dimensionless 

slope,  dimensionless;  also  delta  (impulse) function 

dimensionless  small  quantity,  0  <  e  <<  1 

absolute  viscosity,  kilogram  mas s/meter- second; 
also  coordinate  (y-direction) ,  centimeters 

flow  angle,  radians 

thermal  conductivity.  Watts/centimeter  - 
degree  Kelvin 

wavelength,  centimeters 

Mach  angle  defined  by  arcsin(l/M) 

inverse  radius  of  curvature,  1/centimeter 

coordinate  (x-direction) ,  centimeters;  also 
transonic  similarity  parameter,  dimensionless 

density,  kilogram  mass/cubic  meter 

summation 

standard  deviation  of  Gaussian  beam,  centimeters 

thickness  ratio,  dimensionless;  also,  vibrational 
relaxation  time,  seconds 

viscous  dissipation  function,  Watts/cubic  centimeter 
velocity  potential,  dimensionless 
viscous  force  vector.  Newtons/cubic  meter 
angular  velocity,  radians/second 


Superscripts 

(  ) '  prime,  dimensionless  quantity 

(  )  bar,  normalized  quantity;  also  average  value  of  a 

quantity 

(~)  tilde,  perturbation  quantity 


Subscripts 


a,b 


i 

j 

L 

0 


w 

00 


indicate  values  of  a  quantity  in  front  of  and 
behind  a  shock  wave 

x-direction  or  streamline  index  -  first  subscript 
(i  » 1  denotes  freestream  conditions) 

y-direction  or  normal  index  -  second  subscript  - 
(j=l  denotes  centerline  conditions) 

linear  value 

total  or  stagnation  condition;  also,  peak  value  of 
an  intensity  distribution 

wall  or  surface  value 

freestream  condition 


1,2,  first,  second,  etc,  perturbation  quantities 


Miscellaneous 

V  Del  operator,  1/centimeters 

Laplace  operator,  1/square  centimeters 
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I. 


INTRODUCTION 


The  optical  quality  of  a  propagating  laser  beam  is  de¬ 
graded  when  there  are  density  gradients  in  the  medium.  A 
laser  beam  propagating  through  an  absorbing  medium  creates 
local  changes  in  the  density  which  affect  the  index  of  refrac 
tion  of  the  medium.  Therefore,  the  density  gradients  will 
refract  the  light  into  regions  of  higher  index  of  refraction 
(higher  density)  causing  defocusing  of  the  beam.  This  pro¬ 
cess,  which  is  known  as  thermal  blooming,  has  been  studied 
by  many  authors  [Refs.  1,  2,  and  3]  to  determine  its  effect 
on  a  laser  beam  propagating  through  various  media. 

When  a  laser  beam  is  slued,  a  complicated  coupled  pro¬ 
cess  of  thermal  expansion  and  fluid  flow  (relative  winds) 
takes  place.  Depending  on  the  rate  of  sluing,  the  beam  will 
experience  subsonic,  transonic,  and  supersonic  winds  at  vari¬ 
ous  radial  distances  from  the  laser  source.  The  subsonic 
and  supersonic  thermal  blooming  problems  have  been 
solved  [Refs.  4-8].  A  computer  program  (DIDER)  has  been 
developed  and  experimentally  verified  [Ref.  9]  for  beam 
interaction  with  supersonic  flows  internal  to  a  Gas  Dynamic 
Laser. 

A  solution  for  the  transonic  regime,  and,  in  particular, 
the  case  when  the  freestream  Mach  number  is  precisely  unity, 
has  been  of  considerable  interest  because  the  addition  of 
heat  can  lead  to  extremely  strong  density  gradients.  As 


with  the  theoretical  investigation  of  transonic  flows  over 
airfoils  and  wedges,  the  linear  small  perturbation  flow 
equations  valid  for  subsonic  flow  (elliptic  equations)  and 
supersonic  flow  (hyperbolic  equations)  are  no  longer  complete 
in  the  transonic  regime. 

Many  authors  [Refs.  10  -  15]  have  studied  various  methods 
of  solving  the  transonic  flow  equations  without  heat  addition. 
The  assumptions  based  on  known  transonic  behavior  about  air¬ 
foils  and  wedges  are  used  to  simplify  the  governing  equations 
and,  further,  to  obtain  approximate  numerical  solutions. 

Often,  the  assumptions  greatly  restrict  the  types  of  transonic 
flows  that  can  be  solved.  For  example,  in  the  parabolic 
method  [Ref.  16],  the  transonic  flow  equation  is  transformed 
into  a  parabolic  differential  equation.  Since  the  solutions 
to  parabolic  differential  equations  are  well  established  in 
the  literature,  this  method  may  be  used  when  the  acceleration 
or  deceleration  is  approximately  constant  in  the  region  of 
interest  (i.e.  across  the  cord  of  many  airfoils) ;  however, 
this  method  is  not  valid  in  the  region  of  the  stagnation 
point  or  where  the  flow  acceleration  changes  sign. 

With  heat  addition,  the  ;basic  flow  equations  for  transonic 
flow  become  more  complicated.  Zierep  [Ref.  17]  has  solved 
the  case  when  the  Mach  number  is  precisely  unity;  but,  in 
solving  the  equation,  he  has  reduced  the  problem  to  that  of 
one  dimension.  Others  [Refs.  18  and  19]  have  solved  the 
nonsteady  transonic  equation  with  heat  addition;  however, 
the  results  are  valid  only  for  a  short  time  after  beam  turn 


on  (i.e.  no  steady  state  solution) .  Even  for  the  short  time 
solutions,  extrapolation  techniques  between  the  property 
values  at  the  critical  points  on  the  subsonic  side  and  the 
supersonic  side  are  employed  to  obtain  the  values  across 
the  transonic  regime.  Experimental  investigations  are  being 
conducted  [Ref.  20]  to  determine  the  influence  of  transonic 
sluing  on  thermal  blooming.  These  experiments  indicate  that 
a  steady  two-dimensional  solution  is  possible.  Further 
experimentation  will  provide  better  guidelines  for  a 
theoretical  analysis  of  the  problem. 

In  this  thesis ,  transonic  thermal  blooming  in  steady 
two-dimensional  invisid  flow  is  formulated  and  solved  for  a 
Gaussian  heat  release  distribution.  In  particular,  it  examines 
the  problem  when  the  freestream  Mach  number  is  precisely 
unity.  This  is  accomplished  by  using  a  method  developed  by 
Broadbent  [Refs.  21  -  25]  in  which  the  streamlines  through¬ 
out  the  flow  are  adjusted  in  such  a  manner  that  the  heat 
addition  necessary  to  achieve  these  streamlines  equals  the 
required  heat  distribution  and  boundary  conditions.  The 
method  is  exact  since  specifying  the  streamlines  reduces 
the  nonlinear  coupled  momentum  equations  in  natural  (stream¬ 
line)  coordinates  to  simple  algebraic  linear  difference  equa¬ 
tions  that  can  be  "marched"  through  the  entire  flow  field. 
Boundary  conditions  are  the  freestream  properties  upstream 
of  the  heat  release  region  and  those  on  a  bounding  stream- 
tube  sufficiently  removed  from  the  heat  release  region  that 
the  flow  is  considered  isentropic.  The  boundary  conditions 


on  the  bounding  streamtube  are  calculated  by  methods  similar 
to  those  employed  in  determining  the  velocity  and  pressure 
distributions  over  airfoils  of  known  shape.  Since  stagna¬ 
tion  points  do  not  exist  on  the  bounding  streamtube  flow, 
a  modification  to  the  method  presented  by  Sprieter  and 
Alksne  [Ref.  10]  for  airfoils  and  wedges  was  developed. 

In  Section  II,  the  theory  used  in  analyzing  the  transonic 
thermal  blooming  problem  is  presented  and  includes  a  des¬ 
cription  of  the  assumptions  made  during  the  mathematical 
analysis,  a  summary  of  the  pertinent  equations  derived  in 
the  Appendices,  and  a  discussion  of  the  method  of  solution. 

In  Section  III,  the  numerical  results  for  the  Mach  1.0  flow 
with  a  Gaussian  heat  release  distribution  are  presented. 
Lastly,  Section  IV  gives  a  brief  discussion  of  the  conclu¬ 
sions  and  of  the  other  methods  that  show  promise  in  solving 
the  transonic  thermal  blooming  problem. 


II.  THEORY 


A  laser  beam  that  is  slued  through  the  atmosphere  at  a 
constant  rate  of  ft  radians  per  second  will  experience  regions 
of  subsonic,  transonic,  and  supersonic  flow  in  planes  per¬ 
pendicular  to  its  axis  at  various  radial  positions  along  the 
beam.  Fig.  1.  The  region  to  be  investigated  in  this  paper 
is  the  transonic  region  in  which  mixed  subsonic  and  super¬ 
sonic  flow  exists,  as  in  contrast  to  purely  subsonic  flow  or 
supersonic  flow,  and,  in  particular,  the  case  when  the  Mach 
number  is  precisely  unity?  that  is,  at  a  radial  position 
given  by  Eq.  (1) . 


_/ym; 

rM„  =  1.0  “  r*  _  JT 


(X) 


For  this  analysis,  the  laser  beam  is  initially  of  a  given 

axial  symmetric  Gaussian  intensity  sluing  into  an  isotropic, 

.  -3 

quiescent  medium  of  constant  density  p (kgm  m  ) .  Viscosity, 

and  thermal  conduction  are  neglected.  Justification  for 

this  assumption  is  given  in  Appendix  A.  The  heating  effect 

of  the  laser  beam  on  the  medium  is  approximated  by  a  molecular 

relaxation  process  which  is  rapid  enough  that  the  absorbed 

energy  is  instantaneously  transformed  into  heat.  The  energy 

distribution  is  of  the  form  given  in  Eq.  (2) . 


,  2  2. 

Q  =  la  =  IQ  a  exp  C - — — ) 

2a 


(2) 


18 


Subsonic  side 
of  transonic 


Source 


High 

supersonic 


Supersonic 


Supersonic 

side  of 
transonic 


Subsonic^* 


Streamiines 
Wake  of  beam 


Compression 

wave 


FIGURE  1.  Flow  Regimes  Along  a  Sluing  Laser  Beam 
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where  Q  (Watts  cm  )  is  the  time  rate  of  heat  absorbed  per 

-2  ... 
volume;  IQ  (Watts  m  )  is  the  beam  maximum  intensity; 

a  (cm  is  the  absorption  coefficient;  and  the  exponential 
of  the  Gaussian  with  standard  deviation  a  (cm),  which  is 
related  to  the  spot  size  of  the  beam  by  w  -  J2  a,  is 
centered  at  the  origin  of  the  x-y  coordinate  system  trans¬ 
verse  to  the  beam  axis  at  r*. 

The  theory  employed  for  this  investigation  separates 
the  steady  transonic  flow  with  heat  addition  into  two  en¬ 
tirely  different  flow  regimes.  Fig.  2.  The  internal,  heat 
addition  region  can  be  treated  exactly  by  the  Broadbent 
method  [Ref.  18] ,  in  which  the  nonlinear  equations  of  motion 
and  the  energy  equation  are  solved  numerically  in  natural 
coordinates.  Then,  a  mesh  of  streamlines  and  normals  is 
constructed  throughout  the  flow  field.  Fig.  3.,  which  reduces 
the  momentum  equations  to  simple  algebraic  equations  solvable 
for  the  flow  properties.  Finally,  the  mesh  of  streamlines 
is  adjusted  until  the  properties  calculated  from  the  linear¬ 
ized  momentum  equations  satisfy  the  energy  equation  and 
boundary  conditions.  Boundary  conditions  are  obtained  from 
the  external,  isentropic  flow  region  by  approximating  the 
transonic  solution  over  a  bounding  streamtube  to  the  internal 
region. 

The  governing  equations  in  the  internal  region  have  been 
derived  for  flow  with  heat  addition,  but  without  viscosity, 
heat  conduction  or  body  forces.  A  small  perturbation  analysis 


* 


* 


4 


y 


FIGURE  2.  Transonic  Flow  Configuration  with  Heat  Addition 


FIGURE  3.  Natural  Coordinate  System  with  Flow  Mesh 
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for  transonic  flow  (Appendix  A)  shows  that  viscosity,  heat 
conduction,  and  body  forces  can  be  neglected  to  first  order. 
The  small  perturbation,  steady  two-dimensional  transonic 
equation  with  head  addition  is  written  as: 

(1-M^  )<}>X«XI  +  ^y«yi"Moo  (y+1)$xi$xixi  =  CpT^ 

(3) 

where  the  term  on  the  right-hand  side  is  added  for  flows 
with  heat  addition;  q(x',y*)  (Joules  kgm  is  the  heating 
function;  Cp  (Joules  kgm  ^  °K~^)  is  the  heat  capacity; 

(°K)  is  the  freestream  static  temperature;  is  the 
frees tream  Mach  number;  and  <f>  is  the  nondimens ional  velocity 
potential.  Therefore,  in  vector  notation,  the  steady  state 
equations  [Ref.  26]  are: 


V  (pU)  =  0 

(4) 

pU’VTT  +  Vp  =  0 

(5) 

pTJ"  V  (h  +  =  Q  =  la 

(6) 

p  =  pRT 

(7) 

—  -1  -2 
where  U  (m  sec  )  is  the  velocity?  p  (newtons  m  )  is  the 

static  pressure?  h  (Joules  kgm  1)  is  the  enthalpy?  and 

R  (Joules  kgm  1  °K)  is  the  universal  gas  constant.  Using 

a  two-dimensional  natural  coordinate  system  [Ref.  26] , 

shown  in  Fig.  3,  the  governing  equations  were  derived  from 
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the  above  equations.  The  momentum  equations  are  represented 
in  the  following  convenient  numerical  form  (Appendix  B) . 


-4 


uj 


-  u! 


i+l,j+l  i/j+1 


An ' .  ....  . 

+  — \ /J  (p» 


yu 


i+lfj+l  ^i,j+l^  ® 


(8a) 


-P? 


yM 


(U.',  +U.V,  -V'  •  )  =  0 


*i+l , j+1  *i,j+l  ^2^'  i+l, j+l*i+l, j+1  i+l,  j  i+l, j 

(8b) 

where  ah!  .is  the  streamtube  width  at  the  i,jth  station; 

Y  is  the  ratio  of  specific  heats?  v?  .  is  the  inverse  curva- 

J 

ture  defined  as  AnJ  ./R?  and  the  prime  (')  denotes  dimension- 

J-fj 

less  quantities  with  respect  to  the  freestream  values.  The 
quantities  with  a  bar  and  parenthesized  subscripts  represent 
average  values  between  the  indicated  stations  (i.e.  An^  i+l)  j 
is  the  average  streamtube  width  between  station  i, j  and  i+l, j) , 


and  the  properties  with  the  tilde  (~)  denote  perturbation 

values  from  freestream  conditions  (i.e.  u.  .  =  (u.  .  -  u,  .)) . 

r,  j  if  j  J-f  J 

The  natural  coordinates  s  and  n  are  tangent  and  normal  to 
the  streamlines,  respectively.  Eqs.  (8a)  and  (8b)  represent 
two  nonlinear  equations  since  u',  p*.  An’  and  v'  are  dependent 
variables.  However,  specifying  the  streamlines  of  the  flow, 
these  equations  become  linear  equations  solvable  for 

Qi+l,j+l  and  pi+l, j+l  at  P°int  e  in  terms  of  "i+l, j'  Pi+l,j 

at  point  d  and  u|  j+i/  j+i  at  P°int  Fig.  3. 

The  remaining  flow  properties  j+1  and  T|+^  can 

be  calculated  from  the  continuity  equation,  Eq.  (9) ,  and 
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equation  of  state,  Eq.  (10) ,  derived  in  Appendix  B.  How¬ 
ever,  the  new  properties  at  the  point  e  (i+l,j+l)  may  not 
satisfy  the  energy  equation,  Eq.  (11) .  If  it  is  not  satis¬ 
fied,  an  iteration  process  is  performed  by  varying  the 
streamline  distribution  and  recalculating  the  properties 
until  coincidence  exists  throughout  the  entire  flow  field. 
Hence,  the  solution  of  the  problem  is  reduced  to  finding 
the  streamline  distribution  that  satisfies  Eqs.  (8)  through 
(12)  with  the  boundary  conditions  determined  from  the 
external  region. 


=  1 


p i+1, j+1  i+1, j+1  i+1, j+1 
pi+l,j+l  =  pi+l,j+lRTi+l, j+1 
*i+l,j+l  ‘  k,  j+1  +  (ui+l,  j+l 


(9) 

(10) 


(IaAn 1 ) 


(11) 


d'i+D'j  Asi, j+i 


where 


Qco  = 


P^CpT.U+fW-lJM/) 


An 


1#  j 


(y-l)An,  . 

*•>  J 


(12) 


Insight  into  the  behavior  of  the  streamlines  at  the 
sonic  speed  was  obtained  from  the  solution  of  the  linearized 
small  perturbation  equation  with  heat  addition  [Refs.  4- 
8] .  A  computer  program  (BLOOM)  was  devised  (Appendix  G) 
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which  gave  solutions  for  a  flow  very  slightly  below  the  sonic 
speed  (M^  =  0.999),  Figs.  4  and  6,  and  very  slightly  above 
the  sonic  speed  (M^  =  1.001),  Figs.  5  and  7.  Computation 
of  the  nonlinear  term  in  Eq.  (3)  from  the  linear  results 
indicate  that  the  linear  equations  are  valid  to  Mot  =  0.9999 
subsonically  and  =  1.0001  supersonically.  From  the  linear 
solutions,  a  trend  of  the  streamline  shape  is  obtained  in 
the  near  sonic  regime,  even  if  not  at  sonic  conditions. 

Nature  usually  exhibits  a  strong  propensity  to  accomplish 
changes  in  a  rather  smooth  fashion  (even  a  shock  wave  is  a 
smooth  change  in  conditions  if  viewed  microscopically 
enough);  thus,  a  knowledge  of  flow  conditions  in  the  near 
sonic  range  gives  a  good  clue  to  the  sonic  flow  behavior. 

In  Figs.  4  through  7,  the  Gaussian  heat  distribution  is 
centered  at  the  origin  on  the  centerline,  and  the  relative 
airflow  is  from  left  to  right. 

From  the  appropriate  density  calculations  and  plots, 
discussed  later  in  the  thesis,  it  can  be  seen  that  the  wake 
term  becomes  less  significant  compared  with  the  source  term, 
which  is  directly  associated  with  the  horizontal  velocity 
perturbation,  as  sonic  velocity  is  approached  from  either 
direction  (i.e.  M  <  1  and  M  >  1) .  Figs.  4  and  5  show  the 
horizontal  (longitudinal)  velocity  perturbation  at  Mach 
numbers  of  0.999  and  1.001,  respectively.  Figs.  6  and  7 
show  the  vertical  (transverse)  velocity  perturbation  at  these 
same  Mach  numbers.  Note  that,  subsonically,  the  horizontal 
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FIGURE  4.  Velocity  Perturbation  (u1)  in  Flow  Direction  for  =  0.999 


-20  -15  -10  -5  0  5  10  15  20  25 


FIGURE  5.  Velocity  Perturbation  (u1)  in  Flow  Direction  for  = 1.001 


perturbation,  Fig.  4,  decreases  to  a  negative  valley 
slightly  upstream  of  the  laser  beam  center  and  then  changes 
rapidly  to  an  equal  positive  peak  slightly  farther  downstream 
of  the  beam  center.  The  perturbation  contours  are  of  very 
small  curvature  near  the  beam  center;  the  supersonic  per¬ 
turbations,  Fig.  5,  show  a  negative  valley  near  the  beam 
center  with  a  return  to  freestream  conditions  downstream  of 
the  center.  In  the  supersonic  case,  the  contours  are  nearly 
normal  to  the  flow  throughout,  in  keeping  with  the  character 
of  supersonic  flow;  whereas,  in  the  subsonic  case,  the  con¬ 
tours  are  closed,  in  keeping  with  subsonic  flow  behavior. 

The  patterns  of  disturbance,  however,  appear  to  be  consistent 
as  the  sonic  speed  is  crossed. 

The  vertical  or  transverse  velocity  perturbation.  Figs. 

6  and  7,  show  similar  behavior  in  the  two  speed  regimes. 

As  expected,  the  subsonic  flow  perturbation  profiles  are 
closed  while  those  for  supersonic  flow  are  open  to  infinity. 
In  both  cases,  the  perturbations  exhibit  an  odd  symmetry 
about  the  main  flow  centerline.  For  the  subsonic  case,  the 
perturbations  are  greater  than  zero  to  the  left  of  the  rela¬ 
tive  flow  and  less  than  zero  to  the  right.  This  is  a  direct 
indication  of  the  spreading  of  the  streamlines  as  the  heat 
release  region  is  approached  by  the  air  flow.  Downstream 
of  the  laser  beam,  the  streamlines  again  become  parallel 
with  the  undisturbed  flow  direction.  The  supersonic  case 
shows  the  same  patterns  except  that  the  spreading  effect  on 
the  streamlines  is  felt  infinitely  far  in  the  transverse 


direction,  roughly  along  the  characteristic  lines  from  the 
disturbances . 

The  effect  of  spreading  of  the  streamlines  as  the  beam 
is  passed,  which  can  be  constructed  from  these  plots  of 
longitudinal  and  transverse  velocity  perturbations,  is  thus 
seen  to  be  quite  similar  immediately  below  and  immediately 
above  the  sonic  speed.  Fig.  2  is  a  sketch  of  the  stream¬ 
lines.  It  is  expected  therefore  that  the  streamline  behavior 
at  the  sonic  speed  will  not  vary  radically  from  this  pattern. 

Based  on  the  streamline  distribution  in  Fig.  2,  the 
following  relationship  between  area  change  (streamtube  width 
variation)  and  change  in  total  temperature  (heat  input)  was 
chosen,  Eq.  (13)  . 


dA  _  dn 
A  n 


|(y+x)  (1+b) 


(13) 


where  Tq  is  the  local  total  temperature  of  the  flow  and  $ 
is  a  variable.  One  dimensional  influence  coefficients  with 
area  change  and  heat  addition  [Ref.  27],  Eq.  (14),  were 
used  to  determine  the  required  $  at  precisely  Mach  1.0. 
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G(x)  =  M2(l+-j(y-l)M2)  (-2-dln(n) 
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In  order  for  Eq.  (14)  to  remain  finite  and  to  be  defined 
at  Mach  1.0,  G(x)  must  tend  to  zero  in  the  limit  as  the  Mach 
number  approaches  unity.  If  Eq.  (13)  is  substituted  into 
Eq.  (14)  and  the  limit  taken,  it  is  found  that  it  must  be 
zero  at  precisely  Mach  1.0  (Appendix  C) .  Hence,  the  stream- 
tube  area  change  is  related  to  the  heat  input  by  Eq.  (15) . 
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j 


dT' 

f(Y+l)  -5?^ 

°l,j 


(15) 


where  the  maximum  difference  between  n!  .  and  nl  .,  and  T' 

irl  1/3  o.  . 

— 6  '  , ' 
and  T*  is  of  the  order  of  10  and,  hence,  the  frees tream 
o,  • 


values  can  be  used  in  the  denominator.  Using  this  approxi¬ 


mation  for  the  streamtubes,  an  explicit  formulation  of  the 
streamtube  shape,  at  any  position  s!  .  along  a  streamtube, 

1 1  j 

was  obtained  (Appendix  C) . 
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where  for  a  Gaussian  heat  input  at  Mach  1.0  (Appendix  C) , 
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Once  the  streamtube  shape  is  known,  the  slope  (dy/dx) , 
Eq.  (17) ,  and  the  radius  of  curvature  (R) ,  Eq.  (18) ,  can 
be  calculated  at  every  mesh  point  throughout  the  flow. 
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Therefore,  the  proportionality  between  the  streamtube 
shape  and  heat  release  distribution  gives  not  only  an 
explicit  formulation  for  the  streamtubes  and,  hence,  linear¬ 
ization  of  the  nonlinear  momentum  equations  but  also  simul¬ 
taneous  agreement  between  those  properties  calculated  from 
the  momentum  equations  and  the  energy  equation.  The  problem, 
in  the  case  of  precisely  =  1.0  flow,  becomes  one  of 
determining  the  correct  boundary  conditions. 

The  boundary  conditions  that  are  used  are  those  corres¬ 
ponding  to  frees tr earn  conditions  (u,'  .  =  pi  .  =  pi  .  =  Ti  ,  -  0) 

1,J  X,J  J.,j  x ,  j 

which  are  far  enough  upstream  that  the  properties  are  not 
affected  by  the  perturbation  caused  by  the  heat  release 
region,  and  the  pressure  (p!  .*)  and/or  velocity  (u!  .*)  on 

the  j*th  streamtube  which  is  considered  far  enough  removed 
from  the  heat  release  region  that  the  flow  can  be  considered 
isentropic.  Eqs.  (16b)  and  (16c)  can  be  used  to  determine 
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the  streamtube  shape  and  displacement  where  j  =  j*.  However 

a  more  convenient  form  is  derived  in  Appendix  C  for  a 

Gaussian  heat  release  distribution  at  M  =1.0  which  can  be 

00 

written  as: 
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and  the  slope  of  the  bounding  streamtube: 
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where  t*  is  the  maximum  growth  in  the  streamtube. 

Since  the  mass  flow  is  constant  in  a  streamtube,  the 
bounding  streamtube  can  be  considered  as  a  solid  wall  simi¬ 
lar  to  the  surface  of  a  wedge  with  a  cusp  extending  forward 
to  the  freestream  conditions  upstream  at  infinity.  Fig.  2. 
By  hypothesis,  this  flow  is  isentropic  and  can  be  treated 
by  methods  similar  to  those  used  in  analyzing  transonic 
airfoil  and  wedge  flows.  In  particular,  an  approximate 
solution  is  obtained  through  an  iteration  process  on  the 
integral  equation  derived  by  Green's  functions  for  small 
disturbance  transonic  flow  similar  to  that  performed  by 
Sprieter  and  Alksne  [Ref.  10  andH]  for  thin  airfoils  and 
wedges  but  modified  for  the  streamtube  configuration  which 
has  no  stagnation  points .  The  absence  of  stagnation  points 


in  the  case  of  the  flow  around  the  bounding  streamtube 
enhances  the  accuracy  of  the  approximate  numerical  tech¬ 
niques  because  the  stagnation  point  represents  a  singularity 
in  the  problem  and  places  added  constraints  on  the  method  and 
accuracy  of  the  solution  in  the  region  of  the  stagnation 
point.  In  this  method,  the  quadratic,  nonlinear  nature  of 
the  governing  equations  is  retained  which  precludes  shock- 
free  supercritical  flows  in  the  transonic  regime. 

The  integral  equation  for  transonic  flow  is  derived  by 
a  Green’s  function  analysis  on  the  small  perturbation 
transonic  equation  without  heat  addition,  Eq.  (3)  without 
the  right-hand  side,  resulting  in  Eq.  (21) . 


(24) 


36 


and  the  normalized  variables  are  defined  in  Appendix  D. 

Eq.  (21)  is  a  function  only  of  x,  the  normalized  coordinate 
along  the  wall  surface,  and  T,  the  normalized  amplitude  of 
the  wall  slope,  given  by: 


—  -  t  exp  - -L  (25a) 

dx  2crz 


where 


_  M  2  (y+1)  t  Mj(y2-l)XaJZn  o 

t  =  - * - 2 -  (25b) 

B 

The  iteration  technique  (Appendix  D)  is  performed  for 
various  subsonic  freestream  Mach  numbers  starting  at  the 
initial  supercritical  flow  and  proceeding  toward  the  Mach 
1.0  flow.  Each  solution  obtained  by  the  iteration  technique 
gives  a  unique  T,  which  determines  the  freestream  Mach  number 

and  its  corresponding  normalized  velocity  distribution 

—  —  —  —2  /  3 

Vx'0)  along  the  surface,  if  (u^l)  versus  i  is  plotted 

for  each  position  along  the  surface,  it  is  found  that  as 

the  freestream  Mach  number  tends  toward  unity  the  slope 

becomes  constant.  This  corresponds  to  the  phenomenon  of  the 

Mach  number  freeze  [Ref.  26]  wherein  the  local  Mach  number 

is  invariant  with  changes  in  the  freestream  Mach  number 

when  the  latter  is  near  unity  or,  more  precisely. 
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Vincenti  and  Wagoner  [Ref.  28],  Liepmann  and  Bryson  [Refs. 
29  and  30] ,  and  others  have  shown  that  the  corresponding 
approximate  relation  yielded  by  the  small-disturbance 
transonic  theory  is  given  by: 

=  0  (27) 

i . 

where 


e 


(l-M2 
(Mot2(y+1)  t) 


) 

1/3 


(u-1) 

T271 


(28) 


The  freeze  of  Mach  number  extends  over  a  finite  range  near 
Mach  1.0  where  it  is  constant.  Hence,  £,  the  slope  of  the 
(u-1)  versus  7  curve  at  each  location  x,  is  constant  near 
the  sonic  condition,  and  the  local  Mach  number  at  each 
position  can  be  calculated  by  Eq.  (28)  .  With  the  local 
Mach  number  distribution  known,  the  velocity  perturbation 
can  be  calculated  from  Eq.  (22) ,  and  the  other  flow  proper¬ 
ties  can  be  determined  from  isentropic  flow  relations  [Ref. 
27]  . 

The  Mach  1.0  flow  is  completely  specified.  Boundary 
conditions  are  known  upstream  of  the  heat  release  region 
and  on  a  bounding  streamtube.  Calculation  of  the  entire 
flow  field  is  readily  calculated  from  the  simple  algebraic 
equations  obtained  by  Broadbent's  method,  where  the  stream- 
tube  configuration  is  known  precisely  at  Mach  1.0. 


III.  NUMERICAL  RESULTS 


The  computed  flow  properties  are  for  the  Mach  1.0  thermal 
blooming  problem  with  the  characteristics  listed  in  Table  I. 


Gaussian  Beam 

■  '  ■  ■ 

Peak  Intensity  (IQ) 

3.333x10^  Watts/m^ 

Standard  Deviation  (a) 

5  cm. 

Atmospheric  Absorption  (a) 

8.0  x  10  ^  cm  ^ 

Atmospheric  Freestream 
Conditions 

Temperature  (T  ) 

288.0  °K 

Pressure  (p  ) 

5  2 

1.013x10  Newtons/m 

Density  (pj 

1.2246  kgm/m^ 

Velocity  (u  =  a*) 

340.0  m/sec 

Streamtube  Width 
(An1(j=AnJ 

1  cm. 

Streamline  Cell  Length 
(Asi,i> 

1  cm. 

Ratio  of  Specific  Heats  (y) 

1.4 

Specific  Heat  at  Constant 
Pressure  (Cp) 

1005.0  Joules/kgm°K 

Table  I 

Freestream  Flow  Properties  And 
Laser  Beam  Characteristics 
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The  first  step  in  the  procedure  is  to  determine  the 
boundary  conditions  for  the  heat  release  region.  Far  up¬ 
stream  of  the  heat  release  region,  all  the  dimensionless 
perturbation  quantities  are  zero.  The  other  boundary  con¬ 
ditions  are  calculated  by  Sprieter*s  method  [Ref.  10]  for 
the  bounding  streamtube,  the  shape  of  which  depends  upon 
the  heat  release  distribution  and  is  known,  Eq.  (16) . 

Using  Sprieter's  numerical  iteration  procedure  on  the 
transonic  flow  integral  equation,  Eq.  (21) ,  the  normalized 
velocity  distribution  along  the  bounding  streamtube  was 
obtained  for  several  frees tream  Mach  numbers  between  the 
initial  supercritical  flow  and  Mach  1.0.  The  results  are 
graphed  in  Fig.  8.  From  these  solutions,  a  (u-1)  versus 
t2/3  curve.  Fig.  9,  was  plotted  for  several  positions  along 
the  streamtube.  The  results  indicate  that  the  slope  (£) 
is  constant  for  each  position  and,  hence,  the  Mach  number 
can  be  considered  frozen. 

From  the  values  of  £(x),  the  local  Mach  number  variation 
at  a  freestream  Mach  number  of  unity  was  calculated  from 
Eq.  (28)  ,  since  £  is  known  and  constant  for  the  known 
physical  streamtube  shape,  Eq.  (30) ,  from  Eq.  (D6) . 

T  =  — - —  =  1.386  x  10-6  (30) 

J2t?  a 


where 


t  = 


_  ( y— 1 )  Lairg* 


yp  u 
'  c  00  00 


=  1.737  x  10 


-5 


cm, 
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FIGURE  9 .  Mach  Number  Freeze  £  =  Constant  for 
Various  x  Positions 


is  the  maximum  growth  of  the  bounding  streamtube. 

Once  the  local  Mach  number  was  determined,  the  velocity 
perturbation  along  the  bounding  streamtube  was  calculated 
from  Eq.  (22)  and  is  plotted  in  Fig.  10.  The  remaining 
flow  properties  were  calculated  by  the  isentropic  flow 
relations  [Ref.  27] . 

It  has  been  shown  that  the  area  change  of  the  stream- 
tubes  in  the  heat  addition  region  is  directly  proportional 
to  the  change  in  total  temperature;  that  is,  M  0  in 
Eq.  (13) .  The  energy  equation  is  automatically  satisfied 
by  this  selection  for  the  area  variation.  Since  the  heat 
release  distribution  is  known,  all  of  the  streamtubes  are 
specified.  From  the  streamtube  shape,  the  radius  of  curva¬ 
ture  of  the  streamlines  were  calculated  at  every  point 
throughout  the  region,  Eq.  (18) . 

With  the  curvature  known,  the  calculation  of  the  flow 
properties  in  the  heat  region  is  reduced  to  the  solution 
of  two  simple  linear  algebraic  equations,  Eqs.  (8a)  and 
(8b),  the  momentum  equations  in  natural  coordinates.  The 
results  of  this  "marching”  technique  are  shown  in  Fig.  11 
where  the  density  perturbations  in  the  upper  half  of  the 
symmetric  laser  beam  are  presented.  This  method  has  the 
distinct  advantage  that  it  can  be  solved  on  an  HP  9830 
computer . 

Figs.  12  and  13  show  the  linear  supersonic  and  subsonic 
density  perturbation  solutions  calculated  from  the  computer 
program  (BLOOM)  developed  in  Appendix  G.  The  results 
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>ensity  Perturbation  (p‘)  for  =  1.0 


FIGURE  12.  Density  Perturbations  (p1)  for 


0.999 


presented  in  Fig,  10  at  Mach  1.0  are  consistent  with  a 
continuous  transition  through  the  transonic  regime.  Although 
the  solution  calculated  is  for  precisely  Mach  1.0,  it  can 
be  logically  hypothesized  that  the  flow  behavior  through 
the  transonic  regime  is  as  shown  graphically  in  Fig.  14. 

The  linear  supersonic  and  subsonic  solutions,  and 
Sprieter's  high  subsonic  solutions  are  in  excellent  agree¬ 
ment  with  this  hypothesis  as  well  as  that  conjectured  from 
the  hodograph  plane  representation  for  transonic  flow 
(Appendix  E) .  Therefore,  there  is  a  steady,  two-dimensional 
solution  for  the  transonic  thermal  blooming  problem  which 
is  finite  and  does  not  result  in  "catastrophic"  defocus ing 
of  the  laser  beam. 
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IV.  CONCLUSIONS 


It  has  been  demonstrated  that  there  is  a  two-dimensional 
steady  state  solution  to  the  transonic  thermal  blooming 
problem  for  a  sluing  laser  beam  with  a  freestream  Mach  number 
of  precisely  unity.  According  to  the  linearized  solutions 
for  thermal  blooming,  the  density  perturbations  become 
infinite  as  a  Mach,  number  of  unity  is  approached.  Due  to 
nonlinear  effects  the  trend  to  infinity  is  cutoff,  and  the 
finite  value  of  perturbation  quantities  has  been  established. 

In  obtaining  the  solution,  it  was  found  that  at  precisely 
Mach  1.0  that  the  streamtube  area  variation  through  the  heat 
release  region  is  directly  proportional  to  the  change  in 
total  temperature  in  that  region.  With  the  streamtube 
shapes  known,  an  exact  solution  for  the  flow  field  can  be 
obtained  by  Broadbent’s  method  [Ref.  21]  with  the  corres¬ 
ponding  boundary  values  for  Mach  1.0  isentropic  flow  far 
from  the  heat  release  region.  There  are  no  restrictions  on 
the  existence  of  shock  waves  nor  on  the  symmetry  or  asymmetry 
of  the  heat  release  distribution,  provided  the  bounding 
streamtube  in  the  isentropic  region  and  the  boundary 
conditions  on  it  can  be  calculated. 

Furthermore,  in  obtaining  the  solution  at  precisely 
Mach  1.0,  it  was  found  that  shock  waves  exist  throughout 
the  transonic  regime.  Fig.  14.  As  the  sonic  condition  is 
approached  from  subsonic  speeds,  a  shock  wave  forms  at  the 


downstream  1/e  (one  spot  size)  size,  of  the  Gaussian  heat 
release  distribution,  where  the  local  Mach  number  first 
becomes  sonic,  and  rapidly  moves  downstream  to  infinity. 

If  the  sonic  speed  is  approached  supersonically,  a  shock 
is  formed  at  the  beam  center.  To  be  consistent  with  the 
results  in  the  subsonic  portion  of  the  transonic  regime, 
it  has  been  hypothesized  that  the  shock  moves  rapidly  up¬ 
stream  to  infinity.  Finally,  the  analysis  in  Appendix  F 
indicates  that  there  is  a  minimum  laser  sluing  rate  that 
must  be  maintained  to  preclude  significant  phase  distortion 
in  the  beam  at  Mach  1.0. 

The  results  of  the  Mach  1.0  thermal  blooming  problem 

for  a  Gaussian  intensity  distribution  are  in  excellent 

agreement  with  the  approximate  results  obtained  by  Ellinwood 

and  Mirels  [Ref.  19] .  Fig.  15a  shows  these  results  where 

the  maximum  density  perturbation  was  calculated  to  be  of 
2/3 

the  order  of  Q  '  .  Q  is  defined  as: 

(y-1) I  or 

Q  =  - 2 —  ,  (31) 

Ypa 

where  r  is  the  radius  (spot  size)  of  a  circular  beam  of 

uniform  heating  intensity  (IQa) .  This  approximation,  when 

applied  to  the  Gaussian  beam  of  Section  III,  gives  a  maximum 

-4 

density  perturbation  of  the  order  of  10  which  is  consistent 
with  that  obtained  in  this  paper,  and  shown  in  Fig.  15b. 
Although  the  density  perturbations  are  very  small  (i.e. 
measured  in  parts  per  million) ,  significant  optical  defocusing 


Mach  Nunfeer,  M  = 

(a)  Results  of  Ellinwood  and  Mirels  [Ref.  16] 


(b)  Results  of  Linear  Subsonic  and  Supersonic  Solutions  from  DIDER 
[Refs.  4  and  5]  and  of  Broadbents  Method  at  Mach  1.0 


FIGURE  15.  Schematic  Representation  of  the  Steady  State 
Density  Perturbations  for  a  Sluing 
Two-Dimensional  Laser  Beam 


of  the  laser  beam  could  result  (Appendix  F) .  With  an 

increase  in  the  beam  area  due  to  defocusing,  a  large  decrease 

2 

in  the  beam  intensity  I  (Watts/m  )  will  occur. 

Two  other  methods  that  show  promise  in  solving  the  tran¬ 
sonic  thermal  blooming  problem  are  the  finite  difference 
method  developed  by  Murman  and  Cole  [Ref.  12]  for  isentropic 
transonic  flow  which  can  be  extended  to  include  heat  addi¬ 
tion,  and  the  "finite  element  method"  of  solving  partial 
differential  equations  [Ref.  31]  which  is  a  powerful  method 
for  solving  a  wide  range  of  engineering  problems. 

A  Mach  1.0,  steady  state  solution  for  thermal  blooming 
in  a  laser  beam  with  a  given  Gaussian  heat  release  distribution 
has  been  presented  in  this  thesis.  Fig.  11.  Subsequent  inves¬ 
tigations  should  be  performed  to  ascertain  the  Mach  1.0 
solutions  for  various  Gaussian  and  non-Gauss ian  heat  release 
distributions,  beam  intensities  (I) ,  atmospheric  absorption 
coefficients  (a) ,  and  freestreara  flow  properties.  In  this 
manner,  a  correlation  might  be  developed  that  would  be  use¬ 
ful  as  an  engineering  approximation  to  the  density  perturba¬ 
tions  at  various  transonic  Mach  numbers,  vice  a  simple  order 
of  magnitude  linear  extrapolation  as  presented  by  Ellinwood 
and  Mirels  [Ref.  19] .  Although  solutions  are  known  at  Mach 
1.0  and  in  the  linear  subsonic  and  supersonic  regions,  from 
the  computer  program  BLOOM,  it  would  be  premature  to  estab¬ 
lish  a  criterion  that  could  be  used  in  estimating  transonic 
behavior  in  general  before  additional  work  is  completed. 


APPENDIX  A 


DERIVATION  OF  STEADY  TRANSONIC  FLOW 
FOR  A  CONDUCTING  GAS  WITH  HEAT  ADDITION 

Eqs •  (Al)  through  (A4)  govern  steady  compressible  flow 
with  heat  addition  and  constant  properties  [Ref.  32].  These 
equations  include  viscosity  and  heat  conduction.  The 
distribution  of  the  heat  addition  is  given  by  the  function 
QQf(x,y)  and  is  assumed  known.  The  gas  flow  is  assumed 
uniform  at  infinity  (x  =  ±«>)  . 


V-(pU)  =  0  (Al) 

pU*VtI  +  Vp  =  iff  (A2) 

where  ip  =  rjV(V*U)/3  +  rjV^U 


pTU •  Vs  =  -  V*q  (A3) 

where  ,  $  =  ^  (V*U)2  -  n(V2(U-U)  -  (VXU)2  -  2U-V2U) 

2 

and  V-q  =-kV  T  -  QQf(x,y) 


2  2 
Vp  =  pa  Vs/Cp  +  a  Vp 


or 


s  -  s 


°>  =  T/p^"1* 


(A4a) 


exp( 


Cv 


(A4b) 


The  following  dimensionless  variables  are  defined: 


X' 

=  x/£x 

y'  =  y/Z-2 

p* 

8 

Q. 

\ 

Q. 

il 

T '  =  T/T„ 

u-  =  u/u„ 

s 1 

=  s/cv„ 

n'  =  n/n„  =  l 

8 

5^ 

\ 

V 

II 

where  t ^  and  are  characteristic  lengths  in  the  x  and  y 
directions ,  respectively,  and  the  reference  values  are  the 
freestream  conditions.  Substituting  the  dimensionless 
variables  into  Eqs.  (Al)  through  (A4)  and  eliminating  the 
pressure  term  in  the  momentum  equation,  Eq.  (A2) ,  with  the 
equation  of  state,  Eq.  (A4) ,  the  equations  simplify  to  the 
following: 

p»V»*u»  +  u'.V*-p'  =  0  (A5) 

(A6) 


(A7) 

where  the  dimensionless  ratios  are  defined  as  follows: 


m 

Reynolds  number: 

Re  =  p„u»  lx/n- 

Mach  number: 

M  =  U  //yRT 

co  oo'  ' 

% 

Prandtl  number: 

Pr  =  C  tl>. 

*^oo 

Heat  quantity: 

=  QoVe~u„ 

'm 

Scaling  factor: 

L  = 

p*U'  •V*U*  -  - 


p ' T 1 U 1 -V'S*  = 


T'V'p*  p'T'V's'  L  V*  (V*  -tT* )  ^  V,2U' 
m  .  „  2  3  Re  Re 


Y(Y-1)M„ 


j  -  |  ( V '  -tJ* )  2  +  V'2(U'  -U) 

-  (V'XU')2  "  2  (U '  •  V '  2U ' )  | 


i  , 


+  +  Y  (Y-l)M002H00f  (x,y) 
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For  flows  in  which  conduction  and  viscosity  can  be 
neglected,  two  flows  will  be  dynamically  similar  when  y, 
(Mach  number  in  the  absence  of  heat  addition)  and 
are  the  same. 

Retaining  the  viscous,  conduction,  and  heat  addition 
terms,  the  governing  equations,  Eqs.  (A5)  through  (A7),  in 
two-dimensions  become : 


3p* 


f)  *  (  1  +  T.  1  )  +  ti*  ^  P  V  4.  T.-vr*  

p  (¥xr  +  L  3yr)  +  THF  +  Lv  3y 


(A8) 


% 


Hu 


x-direction  momentum  equation: 


,  /  ,  3u*  ,  -  ,  3u\ 
P  (U  SIT-  +  Lv  3yT* 


T '  3p 1  _  p’T'  3s'  4  32u‘ 

M  2  3x'  yM  2  3x'  3Re  3x' 2 


32v'  L2  32u 
+  ~ - 

3y' 


3 Re  3x' 3y 7  '  Re 


(A9) 


y-direction  momentum  equation: 


■  /  .  3v 1  -  •  3v 1 v 

p  (u  IxT  +  Lv  dyr)  = 


energy  equation: 


P.T.(U.  §§;  +  lv  §§;> 


LT*  3p  * 

Lp'T* 
YM  2 

1  00 

Ml  + 
3y 1 

4L2  32v' 
3Re  ay . * 

L  32u’ 

,  +  _L 

32v' 

(A10) 

3Re  3x'3y 

T  +  Re 

3x'2 

Y(Y-l)Mm2 

Re 

[*{■ 

,  3u* ) 2 
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L  3u'av 
ax'ay1 

2  , 


3Y 


5 } *  <£■  * 1 


3x 


3u' .  2 

3yr 


Y  ,32T'  .  T 2  32T\ 

+  RePr  +  L  3I1 


ax' 


y  Cy-1) f(x,y) 


3Y 
(All) 
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Since  the  transonic  flow  region  is  of  interest,  a 
transonic  expansion  procedure  similar  to  that  used  for 
incompressible  flow  around  an  airfoil  of  thickness  ratio  t 
[Ref.  33]  will  be  used.  In  the  following  analysis,  the 
heat  addition  will  cause  the  perturbations  from  the  free- 
stream  conditions.  Introducing  the  following  perturbation 
quantities : 


u*  = 


P  = 


+  e4/3  u£  ; 

V*  =  ev| 

+  e4/3  p<  ; 

2 

T*  —  1  +  T*  H 

►  e4/3 

s'  =  1  +  e4/^3  §2  ?  1/Re  -  0(e) 


H^-CMe4/3)  ; 


l1  -  1 


£2  -  e'1/3 


L  =  e 


1/3 


the  first  and  second  order  equations  at  Mm  =  1.0  become: 


First  order  equations: 


3u*  3p,* 

aF-  +  ITx*" 


0 


(A12) 


(Y-l) 


0 


(A13) 


(A14) 


Integrating  Eq.  (A12)  gives  u^  =  -p^  which,  when  substituted 
into  Eq.  (A13) ,  illustrates  that  for  first  order  theory  the 


flow  is  irrotational. 


3v£  9u£ 

3 x '  "  3y  *  ~  0 


(A15) 


Finally,  Eq.  (A14)  gives  the  temperature  distribution  as: 


’i  =  (Y-l)  =  -(Y-l) 


(A16) 


Second  order  equations: 


3u« 

ax7" 


du\  3v^  3pl 

+  Pi  83F  +  W-  +  3^  +  ui  15x*  =  0 


3pi 


(A17) 


8u- 

■5aeT 


3u.j  35,'  n  35;  3pi 

+  ui  Tsr  +  pi  sir  =  “  rr  (sir  +  Ti  ir* 


M. 


YM„ 


3i' 

3x' 


4e 


!/3  32u' 


3x 


,2 


(A18) 


3v’  3v' 

P1  Sx5"  +  U1  '5xr 


M. 


_  3p-J  3pl 

<Ti  W*  + 


_ 1_  3s2  1/3  a2vl 

YM„2  3x' 2 


3S2  . 
3x 


.1/3  9  T1 


=  ye-'  '  - t  +  Y(Y-l)  f (x,y) 


3x 


,2 


(A19) 


(A20) 


Combining  Eqs.  (A17)  through  (A19)  and  using  Eqs.  (A12)  and 
(A16) ,  the  transonic  perturbation  equation  is  obtained. 


3u'  3v* 

(Y+D  ^  ^  - 


46^3  ,  _1/3  32n 


3x 


7? 


3x 


,2 


(y-l) f (x,y) 
(A21) 
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The  first  term  on  the  right-hand  side  of  Eq.  (A21) 

represents  the  effect  of  viscosity,  the  second  term  the 

conduction,  and  the  last  the  heat  addition.  As  can  be  seen, 

1/3 

both  the  viscous  and  conduction  terms  contain  a  e  term 

when  the  Reynolds  number  is  of  the  order  of  e  and,  hence, 

are  small  and  can  be  neglected  under  most  circumstances. 

-2/3 

If,  however,  the  Reynolds  number  is  of  the  order  of  e 
or  smaller,  both  terms  become  important  and  must  be  retained 
in  the  perturbation  equation.  Since  to  first  order  the  flow 
is  irrotational,  the  velocity  potentials  u^  =  3<j>/3x*  and 
v|  =  3<j>/3y'  can  be  introduced  into  Eq.  (A21)  giving: 


~(Y+1)  V+x'*'  ++yV 


(3y+1)  e1/3  , 

3  ?x'x'x' 

+  (y-1)  f(x,,e“1/3y,) 


(A22) 


A  similar  analysis  can  be  carried  out  for  a  freestream 

Mach  number  close  to  unity.  In  this  case,  the  freestream 

2/3 

Mach  number  can  be  approximated  by  Mw  =  1  +  Ke  '  and  the 
transonic  expansion  gives: 


<K  +  (y+1)  "  <t>y,y'  "  3'  *  ^x'x'x’ 


-  (y-1)  f (x\e“1/3y>)  . 


(A23) 


APPENDIX  B 


-i 

BRIEF  DEVELOPMENT  OF  BROADBENT 1 S  METHOD  OF 
r  SOLVING  NONLINEAR  EQUATIONS'  OF  MOTION  WITH 

HEAT  ADDITION 

Broadbent's  method  [Ref.  24]  is  an  inverse  solution  in 
which  the  streamlines,  normals,  and  boundary  conditions 
along  at  least  one  streamline  and  normal  are  specified  and 
the  resulting  flow  field  calculated. 

The  governing  equations  for  this  derivation  are  as 
follows  s 


V* (pU)  =  0 

(Bl) 

pU'Vtf  +  Vp  =  0 

(B2) 

* 

pU*V(h  +  iu-O)  =  Q  =  la 

(B3) 

p  =  pRT 

(B4) 

In  natural  coordinates  [Ref.  26],  these  equations 
become: 


( puAn)  =  0 

(B5) 

S  *  S  - « 

(B6) 

» 4  *  is  ■  “ 

(B7) 

• 

J- 

0“  w  (h  *  5"2’  ■ 

(B8) 
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A  mesh  is  constructed  over  the  flow  field  using  stream¬ 
lines  and  normals  which  are  chosen  in  advance  as  in  Fig.  3. 
The  cells  specify  the  streamtube  boundaries  arid  the  mesh 
lines  the  streamline  slope  and  curvature  throughout  the 
flow  field.  With  the  mesh  specified,  Eqs.  (B5)  through  (B7) 
are  nondimensionalized  using  uniform  freestream  flow 
properties  . ,  p,  . ,  p,  . ,  and  T.  .),  a  constant  uniform 

streamtube  width  upstream  of  the  perturbations  (An£  and 
constant  Cp. 


pi,jui,jAni,j  =  1 


(B9) 
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(BlOa) 
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T  ,;,-T!  +  i  M2(y+1)  (u!2  4j_,  -  u]  ^  ^ )  = 


i+l,j+l  i/j+1  2 


-u*2 

li+l,  j+1  i/j+1' 


(IaAn  * ) . +1)  • 

-  vi/i^i^/J  As*  =  o'  (Bll) 

Si/j+l  gi/j+l  {  1} 


i/3 


RTifj 


(B12) 


where  is  a  freestream  heat  quantity  defined  in  Eq. 

j+1  is  a  dimensionless  heating  term?  and  ^  is  a 
dimensionless  inverse  radius  of  curvature. 


(12)  ; 


The  slope  (dy/dx)  and  the  radius  of  curvature  ( R)  of 

the  streamlines  at  a  point  b  (i+l,j)  are  shown  geometrically 

in  Figs.  3.  Due  to  the  symmetry  of  the  problem,  the  slope 

and  the  curvature  of  the  centerline  streamline  are  zero 

(i.e.  for  all  stations  i  >  1  and  j  =  1) .  Assuming  that 

Ay.  .  =  An.  .  and  Ax.  .  =  As.  .  (i.e.  to  first  order  the 
1,3  1,3  1,3  1,3 

mesh  is  rectangular) ,  and  that  the  streamline  shape  is 
known  explicitly,  the  slope  in  finite  difference  form  is 
expressed  as: 
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6i  i 

J-  r  J 


=  ni+l, j~1/2Ani+l, j 
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ni, j  +  1/2  Ani,i 


(Bl3a) 


or: 
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where 


(Bl4a) 


^,1 


+  ^lAni,!+l  ^i,1  +  ^iA^i,1+1+ 7 


i  >  1,  (Bl4b) 

j  1  2 


The  inverse  of  the  radius  of  curvature  (v)  is  defined  in 
Eg.  (B15)  as: 
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d2y/dx2 

(1  +  (dy/dx)  *)  iri 


(B15) 


By  nondimens ionalizing  and  substituting  the  slope  obtained 
in  Eqs •  (Bl3a)  or  (Bl3b)  into  Eg.  (B15)  ,  the  following 


finite  difference  equation  is  obtained: 
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(B16) 


Because  the  magnitude  of  the  slope  6^  j  is  never  greater 
than  1<T6,  it  can  be  neglected  in  the  denominator  simplifying 
the  expression  for  the  curvature.  Substituting  Eq.  (B13b) 
into  Eq.  (B16)  and  retaining  only  terms  of  first  order,  the 
curvature  becomes 
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If  the  mesh  is  not  specified,  Eqs.  (B9)  through  (B12) 
represent  five  equations  in  five  unknowns  and  are  nonlinear 
since  An*,  p1  and  u'  are  dependent  functions.  When  the 
streamtube  shapes  are  specified.  An'  and  the  curvature  v' 
are  known  quantities  and  the  equations  become  linear  in  p* 
and  u*.  Therefore,  Eqs.  (BlOa)  and  (BlOb)  represent  two 
algebraic  finite  difference  equations  solvable  for  p*  and  u 


63 


throughout  the  flow  field,  provided  suitable  boundary  condi¬ 
tions  are  known.  For  example,  if  the  values  of  p*  and  u' 
are  known  at  points  a  (i, j) ,  b  (i+l,j>  and  d  (i,j+l)  in  Fig. 

3,  the  pressure  and  velocity  perturbations  can  be  calculated 
algebraically  at  point  e  (i+l,j+l).  This  procedure  is  con¬ 
tinued  (“marched")  throughout  the  mesh  until  all  flow  pro¬ 
perties  have  been  determined. 

The  remaining  flow  properties,  density,  and  temperature 
are  then  calculated  from  Eqs.  (B9)  and  (B12) ,  the  continuity 
equation  and  equation  of  state.  Once  this  has  been  accom¬ 
plished,  the  energy  equation,  Eq.  (Bll) ,  is  checked  to  see 
if  the  properties  calculated  from  the  specified  streamline 
shapes  agree  with  the  desired  heat  distribution.  If  it  does 
not  coincide,  an  iteration  procedure  is  performed  in  which  the 
streamline  shapes  are  varied  until  agreement  is  obtained. 

Boundary  conditions  are  needed  along  at  least  one  stream¬ 
line  (usually  a  wall)  and  one  normal  (usually  frees tream 
conditions) .  In  the  transonic  thermal  blooming  problem,  the 
boundary  values  are  chosen  as  the  unperturbed  freestream 
properties  upstream  of  the  heat  release  region  and  the 
velocity  or  pressure  on  a  bounding  streamtube  which  is  far 
enough  away  from  the  heat  release  region  that  it  can  be 
considered  isentropic.  Fig.  2.  With  the  pressure  and/or 
velocity  known  on  two  of  the  boundaries,  the  solution  for 
all  the  flow  properties  can  be  calculated  as  mentioned  pre¬ 
viously.  If  a  coundary  condition  is  known  at  other  points 


in  the  flow  field,  the  resulting  properties  from  the 
"marching”  procedure  must  agree  with  it  or  the  specified 
streamtube  configuration  must  be  altered  until  coincidence 
between  the  boundary  conditions  and  heat  distribution  is 
obtained  throughout  the  flow  field. 
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APPENDIX  C  • 


DERIVATION  FOR  THE  RELATION  BETWEEN  AREA 
VARIATION  AND  HEAT  ADDITION  AT  PRECISELY 
MACH  1.0  FLOW 


1.  AREA  CHANGE  WITH  HEAT  ADDITION  AT  THE  SONIC  POINT 

From  the  influence  coefficients  for  constant  specific 
heat  and  molecular  weight  [Ref.  27],  the  change  in  Mach 
number  with  heat  addition  and  area  variation  is  given  by: 


dM  G(x) 
dx  ,  ,2 


(Cla) 


where 


21  2  2 
G(x)  =  MZ  (l-hj(y-l)  M  j  [a+V*T)—££'~ 


2dx  1 '  (Clb) 


d(lnTQ)/dx  and  d(lnA)/dx  are  assumed  to  be  functions  of  x, 
the  streamwise  coordinate,  and  dx  is  always  positive.  Based 
on  the  relation  between  streamline  shape  and  total  temperature 
change  obtained  from  the  linear  subsonic  flow  and  supersonic 
flow,  Figs.  4  through  7,  the  following  relationship  between 
area  variation  and  total  temperature  change  (heat  addition) 
is  chosen: 


j  / 1 _ * \  i  d( InT  ) 

=  i(Y+l)(l+6(x)) 
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$  is  a  variable  providing  an  additional  degree  of  freedom 
in  determining  the  relationship  at  Mach  1.0.  Substituting 
Eq.  (C2)  into  Eq.  (Clb) ,  G(x)  simplifies  to: 

99  d (In  T  ) 

G(x)  =  M^a+i-U-Dtn  ((M^-l)Y  -  (Y+1)B)  .  (C3) 

Since  the  denominator  of  Eq.  (Cla)  tends  to  zero  in 
the  limit  as  the  Mach  number  approaches  unity,  the  numerator 
must  go  to  zero  to  permit  a  continuous  passage  from  subsonic 
to  supersonic  speeds  (or  vice  versa) .  The  condition,  when 
G(x)  passes  through  zero  simultaneously  with  the  Mach  number 
becoming  unity,  is  designated  the  critical  point  G(x)  and, 
from  Eq.  (C3) ,  is  equal  to: 


lim  G(x) 
M+l 


g*(x)  =  -f(Y+ires 


d(ln  Tq) 


CC4) 


Hence,  for  heat  addition  -d(ln  TQ)/dx  >  0,  3  must  equal 
zero  or,  equivalently,  the  area  change  is  directly  propor¬ 
tional  to  the  heat  addition  (i.e.  change  in  total  tempera¬ 
ture)  for  steady  flow  from  Mach  1.0,  Eq.  (C5)  is: 


A)  —  i(y+l) 
dx  2Kr  ±} 


d(ln  Tq) 

33E - 


(C5) 


2 

At  the  critical  point,  dM  /dx  is  of  the  indeterminate 
form  0/0.  Therefore,  the  limiting  value  of  dM  /dx  at  the 
critical  point  is  found  by  applying  L' Hospital's  Rule  to  the 
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right-hand  side  of  Eq.  (Cla) 


lia  'ad*  $£>*  =  ias^. . 

M+l  dx  dx  -  (dM2/dx)  * 


(C6a) 


or ,  equi va len t ly , 


(C6b) 


2 

Eq.  (C6b)  suggests  that  (dM  /dx) *  can  be  positive  or  negative 
and  the  flow  can  become  either  subsonic  or  supersonic  from 
the  sonic  point.  The  requirements  on  3  follow  from  Eq.  (C6b) 
where  G(x)  is  given  by  Eq.  (C3) . 


(— )  * 
dx 


=  i(Y+i) 


-  ( Y+l)  (^) 


d(ln  Tq) 


Thus,  at  the  critical  point,  the  flow  may  be  continuous 
through  the  sonic  speed  provided  (dG/dx) *  is  negative  and 
is  not  zero.  For  heat  addition,  Eq.  (C7)  can  be  solved  for 
the  critical  value  of  (d$/dx)*,  Eq.  (C8) : 


and  is  graphed  in  Fig.  1C. 


k 


FIGURE  1C.  3  Variation  versus  Mach  Number  Variation 

for  Sonic  Flow 


From  the  sonic  flow  condition,  6  must  be  positive  for 
supersonic  flow  and  can  be  either  positive  or  negative 
for  subsonic  flow.  The  change  in  3  with  Mach  number  squared 
must  be  greater  than  y/(y+l)  for  either  flow  condition, 
supersonic  or  subsonic. 

In  conclusion,  for  sonic  flow  the  area  variation  is 
directly  proportional  to  the  change  in  total  temperature, 

Eq.  (C5) ,  where  the  change  of  3  from  the  sonic  point  must 
satisfy  Eq.  (C8)  . 


69 


2.  STREAMLINE  CONFIGURATION  FOR  MACH  1.0  FLOW  WITH  A 
GAUSSIAN  HEAT  INPUT  DISTRIBUTION 

In  the  first  part  of  this  Appendix  the  streamtube  area 
change  was  shown  to  be  directly  proportional  to  the  change 
in  total  temperature  at  precisely  M^  =  1.0?  see  Eq.  (C5) . 

To  derive  an  expression  for  the  change  in  total  temperature , 
a  slightly  different  form  of  the  energy  equation  is  used, 

Eq.  (C9)  , 


p"  h  (CpV  =  Q  = Ia 


(C9) 


in  terms  of  the  total  temperature  TQ.  Assuming  a  constant 
Cp  and  a  Gaussian  heat  intensity  distribution  in  the  form 
given  in  Eq.  (2),  Eq.  (C9)  becomes: 
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where  n!  .  represents  the  total  displacement  of  the  jth 
1/3 

streamtube  from  the  centerline  ( j^l)  ,  and  Q^,  defined  by 
Eq.  (12),  is  evaluated  at  M^  =1.0. 

Equation  (CIO)  is  now  substituted  into  Eq.  (C9)  to  give 
an  expression  for  a  change  in  total  streamtube  growth  over 
a  streamline  distance  As|  at  Mw  =  1.0. 
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or  more  conveniently  for  computational  purposes: 
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#  . 


for  streamtubes  off  of  the  centerline  ( j  >  1) •  Equations 
(C12a)  and  (C12b)  can  be  integrated  to  give  an  expression 
for  the  total  streamtube  perturbation  from  its  initial 
width  n,1  .  as: 
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An  alternate  form  of  Eqs.  (C13a)  and  '(C13b)  which  can  be 
used  in  Eqs.  (Bl3a)  and  (Bl3b)  to  simplify  the  expressions 
for  the  slope  and  the  radius  of  curvature  of  the  streamlines 
is  as  follows: 
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where  it  has  been  assumed  throughout  that  the  initial 

streamtube  widths  are  uniform?  that  is,  An.1  .  =  1  or 

1  j 

Ani  .  =  0  for  all  j . 

J- 1 3 

An  expression  for  a  streamtube,  which  is  far  enough  from 


the  heat  release  region  that  it  can  be  considered  isentropic. 


is  now  derived.  Carrying  out  the  integration  over  y*  in 


Eq,  (C12a)  and  assuming  that  n!  .  >  4a 1  an  expression  for 


the  slope  of  the  bounding  streamtube  is  obtained: 
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Integrating  Eq.  (C15)  with  respect  to  ds|  gives  an  explicit 
equation  for  the  bounding  streamtube  perturbation  from  its 
initial  position  n{  .*  as: 
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APPENDIX  D 


V  INTEGRAL  EQUATION  DERIVATION  FOR  TREATING 

SMALL  DISTURBANCE  TRANSONIC  FLOW  WITH  SHOCKS 
FOR  AN  INFINITE  WALL  WITH  GAUSSIAN  SLOPE 


1.  GENERAL  DERIVATION  OF  INTEGRAL  EQUATION 

In  this  appendix,  an  integral  equation  is  derived  for 
the  small  disturbance  transonic  flow  with  shocks  on  an 
infinite  wall  with  Gaussian  slope.  The  small  perturbation, 
steady,  two-dimensional  transonic  equation  without  heat 
addition  can  be  written  as: 


(1“M  >W+Vy 


=  (1— M  2-M  2 

OO  CO 


^x*  ^  ^x’x*  +  *y 1  0 

(Dl) 


where  the  local  Mach  number  is  approximated  by  the  freestream 

2 

Mach  number  plus  a  perturbation  quantity  (M  (y+l)<b  ,). 
Equation  (Dl)  is  valid  only  in  regions  where  the  necessary 
derivatives  exist  and  are  continuous.  Therefore,  the  integral 
equation  derivation  includes  an  equation  for  the  transition 
through  a  shock?  the  equation  is  the  classical  relation  for 
the  shock  polar  [Ref.  34],  Eq.  (D2) . 
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where  subscripts  a  and  b  refer  to  conditions  ahead  of  and 
behind  the  shock.  If  a  small  perturbation  analysis  is 
carried  out  on  Eq.  (D2)  similar  to  that  performed  in 
obtaining  Eq.  (Dl) ,  the  following  relation  is  found  for 
the  perturbation  velocity  components  across  the  shock  wave: 

(l-M2)  (u;  -  u£)2  +  (v;  -  v£)  2  -  mJ  (Y+l)  (u;  +  u£>  (S;  -  u£)  2/2 

(D3) 

Equation  (D3)  corresponds  to  the  shock  polar  curve  for 
shock  waves  of  small  strength  that  are  inclined  at  any  angle 
between  that  of  normal  shock  waves  and  that  of  the  Mach  lines. 

The  applicable  boundary  conditions  for  the  integral 
equation  analysis  are  that  the  perturbation  velocity  must 
vanish  far  upstream  of  the  disturbance,  Eq.  (D4) : 


36 


X=“°° 
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(D4) 


and  that  the  flow  must  be  tangent  to  the  wall  surface, 
Eq.  (D5) : 
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where  9z'/9x*  is  the  slope  of  the  wall  surface  in  the  x*  (s') 
direction  which  can  be  approximated  by  Eq.  (C15) ;  and  t  is  an 
equivalent  thickness  ratio  for  the  surface  defined  as  the 
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maximum  wall  growth  t  (cm)  in  the  y-direction  divided  by  a 
characteristic  length  /2tF  a  for  1.0. 


t  = 


m  a 


t  = 


(y-1)  I0aTra4 


(D6) 


Also,  it  is  necessary  to  prescribe  that  the  direct  influence 
of  a  disturbance  in  the  supersonic  region  proceeds  only  in 
the  downstream  direction. 

Sprieter  and  Alksne  [Ref.  10]  derive  in  detail  the 
integral  equation  by  a  Green’s  function  analysis  for  transonic 
flow  around  bodies  having  subsonic  freestream  velocities,  and 
a  local  shock  discontinuity.  The  nonlinear  term  in  the 
differential  equation  was  maintained.  This  derivation  can  by 
extended  to  the  streamtube  configuration  of  the  boundary 
condition  for  the  internal  region.  For  the  sake  of  brevity, 
only  the  important  equations  and  pertinent  assumptions  will 
be  given  in  the  following  analysis. 

Since  the  object  of  the  analysis  is  to  determine  the 
pressure  perturbation  and/or  velocity  perturbation  on  the 
bounding  streamtube,  an  equivalent  equation  for  the  velocity 
perturbation  can  be  readily  obtained  by  differentiating 
Eq.  (Dl)  with  respect  to  x';  it  is: 


2  ~  2 
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For  the  Green* s  function  analysis,  it  is  advantageous  to 
normalize  Eq.  (D7)  with  the  following  linear  transformation: 

x  =  x*  ;  y  =  By*  ? 

M2  (y+X)  u'  M2to(Y+l)  u' 

u  =  - - -  ;  v  =  - 5 -  (D8) 

fT 

where  6  =  Yl  -  M^ 

In  this 'way,  Eq.  (D7)  reduces  to  the  following: 


a2u  afu 
3x2  3y2 


where  u  can  further  be  defined  by  the  approximation  to  the 

2  ~  2 

local  Mach  number,  Eq.  (Dl)  ,  namely  M  =  (1  +  (y+l)u),  as: 
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(DIO) 


00  00 

From  Eq.  (DIO) ,  it  is  clear  that  for  a  subsonic  freestream 
Mach  number:  u  <  1  when  the  local  flow  is  subsonic,  u  =  1 
when  it  is  sonic,  and  u  >  1  when  it  is  supersonic. 

The  Green* s  function  analysis  is  a  Neuman  problem  for 
finding  a  function  u  in  a  bounded  region  R,  bounded  by  a 
finite  number  of  closed  smooth  contours  C,  such  that  Eq.  (D9) 
is  satisfied  in  R,  u  and  its  first  partial  derivative  are 
continuous  in  R  (shocks  are  excluded  by  judicious  selection 
of  the  bounding  contour)  and  the  normal  derivative  3u/8y 
is  prescribed  on  the  boundary  of  R,  Fig.  ID. 
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FIGURE  ID.  Integration  Region  for  Green's  Function  Analysis 
on  Flows  Past  Bounding  Streamtubes  at 
Transonic  Speeds 


Green's  theorem  states  that: 


/  (uVG«n-GVu*n)  dC  =  //  (uV2G  -  GV2u)  dR  (Dll) 

C  R 


where  the  directional  derivative  on  the  left  side  are 
0  taken  along  the  normal,  drawn  inward,  to  the  curve  C. 
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The  Green* s  function  G  is  constructed  so  that  it  satisfies 
the  Laplace  equation  V  G  =  6  (£-x,r}-y+yo)  with  3G/3n  =  0  on  . 

G(f,n,x,y-y0)  -  £  ln(li)  =  £  ln(R)  (D12) 

where 

/(x-f)2  +  (y-y0-rf) 2 

R  =  . . . . .  . 

/(X-I)2  +  (y+Yo+fT ) 2 

and  X  and  r\  are  the  variables  of  integration  while  x  and  y 
are  the  coordinates  of  a  point  P.  With  the  Green’s  function 
of  Eq.  (D12) ,  Eq.  (Dll)  becomes: 
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where  the  respective  integration  regions  and  contours  are 
designated  in  Fig.  ID.  Integrating  the  remaining  surface 
integral  by  parts  twice  with  respect  to  If,  applying  the 
Green's  function  boundary  condition  on  C1  (17=0  or  y=yQ)  / 
allowing  the  radius  r  of  the  curve  C2  go  to  infinity,  and 
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decomposing  the  line  integral  over  the  shock  wave  into 
components  parallel  to  the  axes  of  the  coordinate  system 
yields  Eq.  (D14) . 


_  3  .  mmm  1  —2  \  /  “*  1  —2  \  9G 

G  —  (u-~u)  -  (u  -  T  u  )  — 

3f  2  2  31 


3  -  1— 2X  1—2  x  3G 

G  — (u-^u  )  -  (u  -  Tu  )  — 

3?  2  2  31 


G  3n  3G 

.  3?  3?J  a 


r  3u  77  3G 

—  “  u  — 

L  3rT  3nJ 


cos  ( 


cos  (n,£) 


nr?)  1 

n  T\  j 


dn 
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where  uLw  represents  the  linear  wall  velocity  u  when  the 
wall  shape  is  specified. 


(  n  3 u.  —  3G 

J  G  —  -  u  — 

00  3n  3n 


n=o 


(D15) 


Sprieter  and  Alksne  [Ref. 10]  in  their  derivation  delineate 
the  advantages  for  integrating  the  surface  integral  over  R 
by  parts  as  was  done  in  arriving  at  Eq.  (D14) .  The  advantages 
are  important  enough  that  they  are  reiterated  here.  First, 
the  double  integral  of  Eq.  (D13)  shows  a  very  strong  influence 
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of  the  velocities  in  the  region  immediately  surrounding  the 

point  P  since  they  are  multiplied  by  ln(R) .  This  influence 

isllargely  nullified  in  the  double  integral  of  Eq.  (D14) 

because  part  of  the  region  has  a  negative  influence  and  part 

has  a  positive  influence.  The  predominant  influence  in  the 

1  —2 

latter  case  is  furnished  by  the  term  j  u  standing  outside 

the  integral.  Next,  the  contribution  of  distant  regions 

is  diminished  in  importance  in  the  double  integral  of 

Eq.  (D14)  since  their  influence  varies  inversely  with  the 

square  of  the  distance,  rather  than  with  the  logarithm  in 

Eq.  (Dll) .  A  further  advantage  is  that  the  value  of  the 

double  integral  of  Eq.  (D14)  is  continuous  through  a  shock 

wave  rather  than  discontinuous  as  is  the  case  with  Eq.  (D13) 

Lastly,  a  point  of  great  importance  in  the  approximate 

solution  arises  from  the  fact  that  the  integration  by  parts 

1  —2 

provides  extra  terms  (those  containing  -j  u  )  in  the 
integrals  along  the  shock  surface  S  which  combine  with  those 
already  present  in  such  a  manner  that  the  contribution  of 
these  integrals  becomes  very  small  when  the  shock  wave 
approaches  a  normal  wave,  as  is  usually  the  case  at  high 
subsonic  speeds. 

Therefore,  u  must  satisfy  Eq.  (D14)  where  u^  is  given 
by  Eq.  (D15)  and  the  shock  relationship  Eq.  (D16) ,  which  is 
the  normalized  form  of  Eq.  (D3) . 

(*a-"b)2  +  <\-\)2  =  7  (^u^a-V2 


(D16) 


In  the  case  of  a  specified  symmetric  wall  shape,  Eg.  (D14) 
can  be  considerably  simplified  as  it  permits  the  introduction 
of  the  boundary  conditions  which  are  specified  at  the  outset 
of  the  problem,  Eqs.  (D17a)  and  (D17b) ,  into  the  integral 
over  the  wall  surface. 
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where  Z  represents  the  reduced  ordinate  of  the  surface,  and 
7  the  reduced  thickness  ratio  given  by: 

M  2  (Y+l)  T 

7  =  - r -  (D18) 


Substituting  the  relationships  of  Eqs.  (Dl7a)  and  (Dl7b)  into 
Eq.  (D15)  completely  determines  at  the  outset  of  the  analysis 
the  linear  wall  velocity  distribution: 
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2.  SIMPLIFICATION  AND  APPROXIMATE  SOLUTION  OF  THE 

INTEGRAL  EQUATION  FOR  TRANSONIC  FLOW 

To  obtain  a  solution  to  Eq.  (D14) ,  some  further 
assumptions  must  be  made.  They  are:  (a)  all  shock  waves 
lie  in  a  plane  transverse  to  the  beam,  and  (b)  the  shock 
waves  are  normal  (i.e.  normal  to  the  local  flow  direction). 

These  two  assumptions  reduce  Eq.  (D16)  to  an  expression, 

Eq.  (D20) ,  which  permits  an  advantageous  introduction  into 
Eq.  (D14)  to  eliminate  the  integral  over  the  shock  S,  Eq*  (D21) . 
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u(x,y)  =  (x,y)  +  i  u2(x,y)  + 

~  -ooQ  3£ 


ln(-i) 


(D21) 


Equation  (D21)  will  be  used  as  the  basis  for  obtaining 
the  velocity  distribution  on  the  bounding  streamtube. 
Approximate  solutions  of  this  equation  are  obtained  numerically 
by  starting  with  an  assumed  velocity  distribution  and 
iterating  until  convergence  is  obtained.  If  mixed  flow 
exists,  the  initial  velocity  distribution  must  include  a 
proper  discontinuity  complying  with  the  shock  relations  of 
Eq.  (D20) . 
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Due  to  the  double  integral  in  Eq.  (D21) ,  an  iteration 
process  would  be  cumbersome  unless  it  could  be  reduced  to 
a  single  integral  by  introducing  suitable  approximations. 
Oswatitsch  [Refs.  35  and  36]  has  shown  that  approximate 
knowledge  of  the  velocity  distribution  in  the  vicinity  of 
the  surface  can  be  expressed  in  terms  of  the  local  coor¬ 
dinate  y,  the  ordinates  of  the  wall  shape  !Z(x),  and  the 
desired  but  unknown  velocity  distribution  uw(x, 0)  on  the 
surface.  This  permits  the  reduction  of  the  double  integral 
to  a  single  integral. 

Oswatitsch  has  considered  such  an  approximate  relation 

in  which  the  velocity  u(x,y)  starts  from  the  value  uw(x, 0) 

at  the  surface  with  an  initial  rate  of  change  given  by  the 

irrotationality  condition,  Eq.  (D22) ,  and  vanishing  at  large 

—2 

distances  as  1/y  . 


3.u 


3v 


The  resulting  relation  takes  the  following  form: 


(D22) 


u(x,y)  =  - 7  (D23) 

(i  +  (y-y0)/b) 

where  b  is  a  function  of  x  satisfying  the  irrotationality 
condition  at  y  =  yQ.  Solving  for  b  using  the  boundary 
conditions  of  Eqs.  (D17a)  and  (D17b)  yields: 
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2u  (x, 0) 

b(x)  =  -  -£■ - 5-  .  (D24) 

d^Z/dx 

Substituting  Eq.  (D23)  into  the  double  integral  of  Eq.  (D21) , 
integrating  with  respect  to  tT  and  setting  y  =  yQ  results  in 
a  simplified  integral  equation  for  calculating  uw(x,  0). 

For  computational  purposes,  the  integral  equation  is  expressed 
as : 

uw(x,0)  =  uLw(x,0)  +  ~u^(xf0)  -  -jl  (D25) 

where 

“  ^<T,°)  T_x  - 

I  =  /  d£  (D26) 

—00 

and  the  function  E  is: 

E  (■£-=-£)  =  E(X)  =  — i— =-=-[i|Xl  (5  -  10X2  +X4)  - 
b  (I+X2)5  2 

(1  -  10X2  +  5X4)  ln|x|  -i(l+X2)  ( 2 5-71X2-X4— X6 )  ] 
Although  I  is  a  function  of  u  and  is  unknown,  it  is  informative 

W 

to  rewrite  Eq.  (D25)  by  solving  for  uw  in  terms  of  I  and  u^, 
thus : 

uw  =  1  ±Jl  -  (2Uj.w-l)  =  1  ±J I  -  L  ,  (D27) 
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where 


L  -  2uLw  '  1 

a  known  function  since  u_  is  known  once  the  wall  shape  is 

Lw 

specified.  It  should  be  noted  from  Eq.  (D27)  that  the 

discriminant  must  always  be  positive  in  order  to  obtain  real 

values  for  u  ,  thus  I  >  L.  Furthermore,  the  choice  of  the 
w  — 

plus  or  minus  sign  determines  whether  the  local  velocities 
are  subsonic  or  supersonic.  A  change  in  sign  at  a  point 
where  the  radical  is  zero  corresponds  to  a  smooth  transition 
through  the  sonic  speed.  A  change  in  sign  at  a  point  where 
the  radical  is  not  zero  corresponds  to  a  discontinuous  jump 
in  velocity,  namely,  a  normal  shock  from  supersonic  to  sub¬ 
sonic  flow  when  progressing  in  the  flow  direction.  A  dis¬ 
continuity  in  the  reverse  direction  is  inadmissable  since 
it  corresponds  to  an  expansion  shock,  an  impossible  phenomenon 
which  violates  the  second  law  of  thermodynamics. 

3 .  COMPUTATIONAL  TECHNIQUE 

In  this  section,  a  brief  description  of  the  computational 
technique  used  to  solve  Eq.  (D25)  will  be  presented.  Sprieter 
and  Alksne  [Ref.  10]  describes  in  great  detail  how  a  judicious 
selection  of  the  initial  velocity  distribution,  depending 
upon  the  surface  shape  and  flow  configuration  (i.e.  entirely 
subsonic  flow  or  mixed  flow  with  possible  shocks) ,  is  essen¬ 
tial  in  providing  a  rapid  convergence  to  the  solution.  It 
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will  be  assumed  for  the  following  discussion,  therefore, 
that  a  judiciously  selected  velocity  distribution  is  known. 
The  method  of  solution  is  as  follows: 

a.  From  the  known  surface  shape  in  normalized  coor¬ 
dinates  Z*,  calculate  its  first  two  derivatives  with  respect 

to  the  normalized  flow  coordinate  x?  that  is,  dsf/dx  and 
2—  —2 

d  Z/dx  .  Since  the  normalization  depends  upon  the  frees tream 
Mach  number,  which  is  unknown,  the  value  of  7,  Eq.  (D18) , 
will  be  a  parameter  in  the  iteration  technique. 

b.  Calculate  the  linear  wall  velocity  uLw(x, 0)  from 
Eq.  (D19)  .  As  7  does  not  depend  on  IT,  it  can  be  taken 
outside  of  the  integral  of  Eq.  (D19)  giving  an  expression 
directly  proportional  to  the  parameter  7. 

^(7,0)  =  7  WLw(7,0)  (D28) 

where  WLw(x, 0)  is  a  function  of  only  the  position  x  and  is 
fixed  once  the  surface  shape  is  specified.  The  quantity 
L  follows  immediately  from  Eq.  (D27) . 

c.  Assume  an  initial  value  of  7  consistent  with  the 
velocity  distribution  and  calculate  b(x)  from  Eq.  (D24) . 

d.  Since  b(x)  is  now  known  and  an  initial  wall  velocity 
uw(x, 0)  has  been  selected,  Eq.  (D26)  can  be  integrated 
numerically  to  get  I. 

e.  The  I  function  calculated  for  the  initial  value  of 
7  is  graphed.  7  is  then  varied  in  the  L  function  until 
the  two  curves  are  tangent  at,  at  least,  one  point  (i.e. 
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I-L  =  0) .  As  was  stated  earlier,  I  >  L  at  all  points  x 
along  the  surface  for  the  velocity  to  be  real. 

f.  Using  the  L  function  that  satisfies  the  tangency 
requirement,  a  new  approximation  to  the  solution  is  obtained 
from  Eq.  (D27) ,  where  the  sign  chosen  for  the  radical  is 
important. 

g.  The  iteration  procedure  continues  with  the  new 
uw(x,0)  and  7  until  both  quantities  converge  to  a  solution. 

h.  Once  a  solution  is  obtained,  the  freestream  Mach 
number  for  a  given  flow  condition  is  calculated  from 
Eq.  (D18)  and  the  local  Mach  number  distribution  from 
Eq.  (DIO) . 


APPENDIX  E 


HODOGRAPH  TECHNIQUES  FOR  TRANSONIC  FLOW 

The  hodograph  plane  represents  a  coordinate  transfor¬ 
mation  in  which  the  velocity  components  are  the  independent 
variables.  Spatial  coordinates  x  and  y  are  replaced  by 
velocity  components  U  and  that  is,  the  physical  coordi¬ 
nates  are  represented  by  x(G,^)  and  y(u,v)  [Ref.  26].  This 
coordinate  transformation  is  convenient  because  it  enables 
the  simple  presentation  of  data  or  solutions  and,  more 
importantly,  it  linearizes  the  two-dimensional  planar  tran¬ 
sonic  flow  equation  that  is  nonlinear  in  the  physical  plane. 
This  equation  is  known  as  the  Tricomi  or  Tricomi-Euler  Equa¬ 
tion.  Boundary  conditions  are  difficult  to  satisfy  [Ref.  15] 
The  small  perturbation  transonic  equation  without  heat 
addition  can  be  written  as: 


(1-M2)|2| 
oo'  3xr 


(y+1)  U 


where  for  ir rotational  flow 


(El) 
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Using  the  hodograph  transformation,  Eqs.  (El)  and  (E2) 
become : 
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u-m^MI  +  Ml  =  M2(Y+i)u'  mi 

3u*.  3v' 


oo '  ~ , 
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(E3) 


and 


Ml  _  Ml =  o 

3v*  3u* 


(E4) 


Eqs.  (E3)  and  (E4)  are  the  linear  transonic  hodograph 
equations . 

As  is  mentioned  in  the  introduction,  the  linearized 
small  perturbation  equations  for  supersonic  flow  are  hyper¬ 
bolic  and  those  for  subsonic  flow  are  elliptic.  The 
transonic  hodograph  equations  change  from  elliptic  to 
hyperbolic  when  the  coefficient  of  3y,/3v'  changes  sign: 


(l-M2)  -  M2  (v+1)  u'  =  0 

00  CO 


(E5) 


The  critical  velocity,  where  the  change  of  sign  occurs, 
can  be  directly  solved  from  Eq.  (E5) . 


u*  * 
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00 
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00 
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Therefore,  for  u'  <  u'*,  Eq.  (E3)  is  elliptic,  and  when 
u*  >  u'*,  it  is  hyperbolic.  The  features  of  transonic  flow 
are  captured  by  the  transonic  hodograph  equations.  Since 
the  flow  occuring  in  the  external  region  and  on  the  bounding 
streamtube  is  considered  to  be  isentropic  two-dimensional 
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V 
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and  planar,  the  hodograph  method  gives  a  conceptual  idea  of 
what  the  resulting  flow  configuration  should  be  over  an 
infinite  wedge  of  prescribed  shape. 

Guderley  and  Yoshihara  [Ref.  IS]  use  the  flow  angle  0 
and  T)r  defined  in  Eq.  (E7)  ,  as  the  independent  variables. 

n  =  (t+l) 1/3  ( (w-w*)  /w*)  (E7) 

r 2 

where  w  =  ^/u  +  v  .  With  this  selection  of  independent 

variables,  the  streamf unction  in  the  hodograph  plane  is: 

*nn  ~  ri  *ee  =  0  (E8) 

and  the  velocity  potential  function  is: 


♦nn  ‘  n  <*’ee::=  0  (E9) 

The  formulation  is  for  a  flow  with  constant  stagnation 
temperature  and,  hence,  can  be  applied  to  the  solution  for 
the  bounding  streamtube.  For  flow  over  a  wedge,  Guderley 
and  Yoshihara  found  a  steady  solution  for  Ma  precisely  unity. 

Fig.  lEa  shows  the  flow  over  the  bounding  streamtube 
in  the  physical  plane.  Points  A,  B,  C,  D,  and  E  are  shown 
mapped  in  the  hodograph  plane.  Fig.  lEb.  A  limiting  Mach 
line  is  shown  in  the  physical  plane.  Any  characteristic  in 
the  region  bounded  by  the  sonic  line  and  the  limiting  Mach 
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line  can  interact  with  the  sonic  line.  Any  characteristic 
downstream  of  the  limiting  Mach  line  does  not  interact  with 
the  sonic  line,  and,  hence,  the  flow  is  not  affected  by  the 
subsonic  flow  upstream  of  the  sonic  line. 
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APPENDIX  P 


DERIVATION  OF  THE  SLUING  RATES  NECESSARY  TO 
PRECLUDE  SIGNIFICANT  PHASE  DISTORTIONS  IN 
A  LASER  BEAM  AT  MACH  1.0 

When  the  phase  difference  between  two  points  in  a  laser 
beam  is  greater  than  a  quarter  of  a  wavelength,  the  deviation 
of  the  beam  from  its  original  direction  is  sufficient  to 
cause  considerable  distortion.  The  range  of  laser  sluing 
rates  necessary  to  prevent  such  phase  distortions  and,  hence, 
to  preclude  defocusing  of  the  beam  at  Mach  1.0  is  derived. 

Consider  a  laser  beam  that  is  slued  as  shown  in  Fig.  1. 
The  Mach  number  of  the  air  relative  to  the  beam  varies  along 
the  beam  as: 


All  Mach  number  regimes  occur  along  the  beam.  If  L  repre¬ 
sents  the  radial  distance  between  two  stations  on  the  beam 
where  the  Mach  numbers  are  M1  and  M2,  it  can  be  determined 
from  Eq.  (FI)  that  L  is  equal  to: 

(M«-M, ) a 

L  =  — V" •  <«> 

The  number  of  waves  (N)  of  laser  radiation  in  the  region 
of  length  L  is  given  by: 
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N 


(F3) 


L  _  n\|L 
X  c 

where  X  is  the  wavelength  of  the  laser  beam  and  is  equal 

to  nv/c.  v  (radians  sec  is  the  laser  frequency  equal 

13 

to  approximately  2  x  10  Hertz  for  a  CO  2  laser,  and  c  (m 
sec"”1)  is  the  speed  of  light.  Therefore,  the  change  in  the 
number  of  waves  in  the  beam  length  L,  due  to  a  change  in 
the  index  of  refraction,  is  obtained  by  differentiating 
Eq.  (F3)  . 

dN  =  — —  (F4) 

c 

The  change  in  the  refractive  index  is  related  to  the 
density  perturbations  in  the  beam,  a  known  quantity,  through 
the  Gladstone-Dale  Law  [Ref.  37],  a  simplified  version  of 
the  Lorentz-Lorenz  Law  for  gases  with  n  =  1,  Eq.  (F5) : 

n  =  1  +  kp  ,  (F5) 


where  k  is  the  Gladstone-Dale  constant  which  is  approximately 
equal  to  2.25  x  10  *  m^kgm  for  air  with  wavelengths  beyond 
the  near  infrared  (X  >  ly) .  Differentiating  Eq.  (F5)  and 
substituting  it  and  Eq.  (F2)  into  Eq.  (F4) ,  the  desired 
expression  for  the  change  in  the  number  of  waves  is  obtained 
as : 
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It  is  of  interest  to  determine  the  range  of  sluing 
rates  for  which  the  change  in  the  number  of  waves  is  less 
than  a  quarter  of  a  wavelength  at  Mach  1.0  from  Eq.  (F7) : 


Q  > 


4k(M2-M1)  avpt 


(o' 

vpMax 


^Min5 


(F7) 


where  is  the  maximum  positive  density  perturbation 

and  p*.  is  the  minimum  (maximum  negative)  density  perturba- 
Mxn 

tion  in  the  beam  at  Mach  1.0*  points  A  and  B  in  Fig.  11 r 
respectively.  As  an  example f  using  M2  —  1.001/-  M^  —  0.999, 
the  flow  properties  in  Table  I,  and  the  maximum  difference 
in  the  density  perturbations  at  Mach  1.0  from  Fig.  11  of 

__  ^ 

1.5  x  10‘4,  the  sluing  rate  must  be  greater  than  7.5  x  10 


radians/sec. 


APPENDIX  G 


SUBSONIC  AND  SUPERSONIC  THERMAL  BLOOMING 
AND  COMPUTER  PROGRAM  (BLOOM) 

Tsien  and  Beilock  [Ref.  38]  have  obtained  solutions  for 
the  linearized  equations  of  motion  with  heat  addition  for 
subsonic  and  supersonic  flow.  The  solutions  are  for  a  line 
heat  source  q  in  a  uniform,  two-dimensional  planar  flow 
which  can  be  considered  to  be  a  perturbation  on  an  initially 
one-dimensional  flow. 

In  1971  and  1972  studies  were  conducted  at  NPS  by  Fuhs 
[Ref.  39]  on  the  applicatin  of  external  heat  to  the  reduction 
of  drag.  These  studies  required  thorough  analysis  of  heat 
addition  in  both  subsonic  and  supersonic  flows*  The  analysis 
of  the  flow  with  heat  addition  found  two  other  applications. 
One  application.was  external  burning  assisted  projectile.  This 
application  was  studied  by  CDR.  W.  J.  H.  Smithey  [Ref.  40]. 

The  other  application  was  recognized  by  Fuhs  [Ref.  4]?  this 
involves  the  density  inhomogeneity  arising  from  quantum 
inefficiency.  Professors  Biblarz  and  Fuhs  [Ref.  5]  developed 
the  capability  and  extended  the  analysis  to  predict  density 
inhomogeneity  in  a  large  gas  dynamic  laser;  see  the  paper 
by  Fuhs,  Biblarz,  Cawthra  and  Campbell  [Ref.  6]  for  details 
of  the  experiment.  Subsequent  to  the  application  of  linear¬ 
ized  solutions  to  flow  internal  to  a  laser,  it  became  apparent 
that  the  analysis  could  predict  flow  properties  in  thermal 
blooming.  In  a  series  of  papers  and  presentations,  Fuhs, 


Biblarz,  Burden  and  Carey  [Refs.  7,  8  and  41-43]  developed 
the  appropriate  analysis  for  linearized  thermal  blooming. 

1.  SUBSONIC  FLOW  WITH  HEAT  ADDITION 

The  perturbations  to  the  flow  [Ref.  38]  are  given  by: 

(Gl) 

(G2) 

(G3) 

P'  - - - y-V"? - «<y)  I(x)  (G4) 
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p'  =  -  jrcBp—  ZF777] 

(y-l) M  q 


where  the  symbols  are  defined  in  the  list  of  symbols  and 
8  =  /l  -  for  subsonic  Mach  numbers.  The  second  term  in 
Eq.  (G4)  is  the  wake  of  the  line  heat  source  q.  The 
temperature  perturbation  can  be  obtained  from  the  equation 
of  state. 

Biblarz  and  Fuhs  [Ref.  6]  in  applying  Eqs.  (Gl)  through 
(G4)  to  laser  internal  aerodynamics  have  presented  the 
integral  equations  used  to  solve  the  unbounded  heat  release 
region  problem  (i.e.  subsonic  thermal  blooming)  and  are 
listed  below  for  completeness.  Since  the  heat  source  is 
volume  distributed,  instead  of  a  single  infinitesimal  rod, 
an  integral  representation  was  used  to  account  for  all 
contributions  to  the  flow  properties  over  the  volume.  The 
geometry  involved  in  the  calculation  is  shown  in  Fig.  1G. 


The  heat  quantity  h(x,y)  is  related  to  the  laser  output 
intensity  through  the  kinetic  lag  that  may  be  present.  To 
develop  an  equation  incorporating  the  lag  between  laser 
energy  and  heating,  it  is  necessary  to  examine  the  kinetics 
of  relaxation.  An  exponential  function  describes  the 
relaxation  with  a  relaxation  time  t [Refs.  44  and  45].  The 
heat  release  function  is  given  by: 

x  ,  *  \ 

h(x,y>  =  la  /  f(x’,y)  exp[-  tx„~  *  '  ]  dx'  (G9) 

where  Q  =  IQa  f(x,y)  is  the  laser  energy  per  time  per 
volume.  The  volume  can  be  interpreted  as  an  area  times  the 
span  of  the  heat  addition  zone  across  the  planar  flow. 

If  significant  variations  of  beam  intensity  in  the  radial 
(beam)  direction  exist,  the  two-dimensional  heat  addition 
function  can  be  approximated  from  the  three-dimensional 
h(x,y,r)  by  an  average: 

ri+l 

h(x,y)  »  —  I  h(x,y,r)  dr  (G10) 

ri 

where  Ar  =  ri+i  ”  ri  ■  T^e  beam  is  divided  into  segments  Ar 
in  length  and  Eqs.  (G5)  through  (G8)  are  solved  for  each 
segment. 

If  the  heat  release  h(x,y)  is  constant,  the  equations 
can  be  integrated  in  closed  form.  The  resulting  integrals 


(Y-I)lo 

S’  =  2,YpJu:  M(x'y;S'n'g) 


(Gil) 


v'  = 


(Y-l)I  a 

°  N(x,y;g,n,B) 


2ryp  Bu 


(G12) 


(y— 1)  MmI  a 
P'  "  -  ZTra^p-— 


(G13) 


(y-l)MIa 

p’  =  -  - 2  M(x,y;g,n,B) 

2*a„  PP» 


(Y“l)  la 

— 2 — —  I  / <S(y-Ti)  Kx-S)  d£dn 
aoo  UooPoo  C  n 


(G14) 


where 

M(x,y?5,n,3)  = 


i^ln 


(x-g)2  +  B2  (y-n) 2 
B2 (y-n) 2 


+  i*z£L  tan-i 


3 (y-n) 

(x-5) 


n 


U  nl 


N(x,y?£,n,3)  = 


( x~  ) 
232 


In 


(x-g)2  +  B2  (y-n) 2 


_  (*-£>  +  (y-n)  t  -i 

„2  +  6 


B 

'  (x-g) 
3 (y-n) 


^2  -*  ^ 


The  first  term  in  Eqs.  (G8)  and  (G14)  correspond  to  the 
density  perturbations  in  regions  outside  the  wake  and  the 
last  term  represents  the  contribution  of  the  wake  and  can 
be  interpreted  as  the  consequence  of  the  heat  convected 
along  the  streamlines  by  the  flow. 
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2.  SUPERSONIC  FLOW  WITH  HEAT  ADDITION 

The  perturbations  to  the  flow  are  also  formulated  by 
Tsein  and  Beilock  [Ref.  38]  as: 


u* 


P’ 

P* 


-  7 F'g'ii3  six  -  gy> 
2YPMeu„ 

(G15) 

ixiila  6 (x  -  gy) 

SYP^u* 

(G16) 

(Y-Dt^q 

Sa  ft—  6(x  -  Sy) 

OO^J^  CO 

(G17) 

(Y-DM^q 

-  -j — . -  <5(x  -  gy) 

_  (Y~Ha —  6  (y)  I  (x) 

(G18) 

- = -  w  /  3 

2a00  3  POT  ^ooaoo  Poo 


Fuhs  [Refs.  4  and  39]  has  modified  these  equations  into 


an  integral  form  as: 


u 


i 


(Y-l)  fS 
2YP00^u„  1 

0 


h(£,n)  sin  U  dS 


v* 

P* 

P' 


„  (Y-l)  f 

2TPoouoo  o 

(H)mm  rs 

=  2CFp—  J 
0 

(Y-1)Mm  s 

■  — 3 —  J 

2a0/3poo  0 


h(£,n)  sin  y  dS 

h(£,n)  sin  y  dS 

h(€,ri)  sin  y  dS 


(G19) 


(G20) 


(G21) 


-  |  J  h(£,ri)  <S  (y-ry)  I(x-£)  d£dn  (G22) 

Mooaco  Poo  5  n 


where  the  symbology  is  the  same  as  that  used  in  the  subsonic 
equations  except  that  8  =  /M2  -  1  for  supersonic  Mach  numbers. 
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3.  COMPUTER  PROGRAM 

The  FORTRAN  computer  program  BLOOM,  included  at  the  end 
of  this  appendix,  was  written  to  compute  and/or  plot  the  flow 
perturbations  (u*,  v',  p'  and  p*)  in  the  subsonic  and  super- 
sonic'  flow  with  heat  addition.  Using  a  modified  version  of  the 
computer  program  DIDER  [Ref.  5]  as  a  subroutine,  the  subsonic 
perturbation  equations,  Eqs.  (G5)  through  (G9) ,  and  the 
supersonic  perturbation  equations,  Eqs.  (G19)  through  (G22) , 
are  numerically  integrated  to  obtain  the  linear  thermal 
blooming  solutions.  The  program  calculates  the  two- 
dimensional  field  at  any  cross  section  of  the  laser  beam 
where  the  energy  release  distribution  is  known.  Finite 
kinetics  are  incorporated  as  described  by  Eq.  (G9) . 

The  program  (BLOOM)  is  versatile  in  that  it  can  calculate 
the  flow  quantities  for  several  subsonic  and/or  supersonic 
freestream  Mach  numbers,  heat  release  distributions,  relaxation 
rates,  and  freestream  flow  properties  in  a  single  computer  run. 
Since  the  length  of  time  for  each  individual  case  depends 
primarily  on  the  size  of  the  flow  region  and  on  the  mesh  that 
is  superimposed  upon  it,  time  requirements  must  be  determined 
on  an  individual  basis.  For  the  results  presented  in  Figs.  4 
through  7  and  Figs.  12  and  13  from  the  program  BLOOM,  the 
two-dimensional  region  of  interest  was  60  cm.  square  centered 
at  the  origin  of  the  Gaussian  heat  release  distribution  and 
was  descretized  into  a  mesh  of  cells  1  cm.  square  for 
subsonic  flow  (M^  =  0.999)  and  into  rectangular  cells  whose 
size  is  determined  by  the  Mach  angle  for  supersonic  flow 
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(M^  =  1.001).  Supersonically,  the  ratio  of  cell  height  Ay 
to  cell  length  Ax  is  equal  to  tan  y.  The  subsonic  results 
take  approximately  25  minutes  and  140  K  of  memory  (without 
plotting  subroutine  CONTUR)  while  the  supersonic  results 
take  approximately  3  minutes  and  285  K  of  memory  (without 
plotting  routine  CONTUR)  to  compile  and  execute  using  an 
IBM  360  computer. 

The  program  can  also  calculate  constant  heat  release 
distributions  for  single  and  multi-pass  laser  beam  configura¬ 
tions  with  or  without  wall  reflections  and  kinetics.  In  its 
present  form,  the  program  can  not  calculate  supercritical 
(mixed)  flows  across  the  beam. 

A  main  feature  of  the  supersonic  portion  of  the  program 
is  that  it  sets  up  the  characteristic  lines  belonging  to  each 
individual  cell  (i.e.  their  diagonals)  and  adds  up  the 
contributions  for  each  flow  property  from  the  characteristics 
and,  in  the  case  of  the  density  perturbation,  from  the  wake 
at  each  cell.  The  calculations  along  the  characteristics 
are  facilitated  by  the  choice  of  cell  shape  and  size  as 
described  above. 


onnn  onnnn  onooo  or»o  nonnnnnoonoooonnooon 


PROGRAM  BLOOM 


COMMON/ONF/T  JM, BTA, ZMA y XDEL * YDELyZI Oy SGM 
COMMON/TWO/ACH,AWAKfXLAG,ZH,ZLT4it 
COMMON/THREE/GMA  fTA, P AyUy ARHOyTAU 

SPECIFY  IN« JN  CONSISTENT  WITH  MATRIX  DIMENSIONS • 

A,  SUPERSONIC  FLOW  (MACH  NUMBER  >  0 )  MllMaCD  nc 

1.  JN  SHOULD  BE  EQUAL  TO  OR  GREATER  THAN  THE  NUMBER  OF 

ROWS  IN  THE  Y-DIRECTION  (JM).  ^  w  r  .M  ,1 

o  FROM  THF  NUMBER  OF  CELLS  IN  THE  Y-DIRECTION  (JM-lIy 
CALCULATE  THE  CELL  WIDTH  ( YDELI  FROM  THE  KNOWN  TOTAL 
Y-DI  MENS!  ON  (ZH);  THAT  IS:  YDFL  -  ZH/(JM-1).  _ 

B.  DIVIDE  YDEL  BY  THE  TANGENT  OF  THE  MACH  ANGLE  TO  OBTAIN 
THE  CFtL  WIDTH  IN  THE  X-DIRECTION  ( XDEL )  • 

4.  DIVIDE  XDEL  INTO  THF  TOTAL  X-DIMENSION  (ZU  JO  DETER¬ 
MINE  THE  NUMBER  OF  CELLS  IN  THE  X-DIRECTION  UM-l). 
DIMENSION  (IN)  SHOULD  BF  AT  LEAST  IM+1. 


B. 


5. 

fVBS3^?SFLE§8Arf§HT^MSoHVR0^  ROWS  I M  THE  Y-DIRECTION 
JM. 

2.  IN  EQUALS  JN. 


I N = 6  5 
JN=65 


SUBSONIC  WITH  CONSTANT  HEAT  INPUT  -  SET  IFL  AG=l 


I FLAG=0 

DENSITY  PERTURBATIONS: 


TABULATED  ONLY 
GRAPHED  ONLY 
TABULATED  AND  GRAPHED 


-  SET 

-  S^T 

-  SET 


JFLAG=! 

JFLAG=2 

JFLAG=3 


JFLAG=3 

VELOCITY  AND  PRESSURE  PERTURBATIONS: 


TABULATED  ONLY 
GRAPHED  ONLY 
TABULATED  AND  GRAPHED 


-  SET  KFLAG=I 

-  SET  KFLAG=2 

-  SET  KFL AG=B 


KFLAG=B 

PHIXX  AND  PHI YY : 


TABULATED  ONLY  -  SET  LELAG=1 
GRAPHED  ONLY  _  “  SET  LFLAG=2 
TABULATED  AND  GRAPHED  -  SET  LFLAG=3 


oononoooon  nnnoonn 


LFLAG-0 


i  |p*. 


c 


’V 


flltUUNITSSHUSr  gENCONSISTFNT  TO  INSURE  DIMENSIONLESS  OUTPUTS. 
THE  UNITS  USED  FOR  THESE  PARTICULAR  RUNS  WERE  AS  FOLLOWS: 

GMA  -  DIMENSIONLESS  % 

TA  (AMBIENT  TEMPERATURE)  -  DEGR.EES  KELV IN 
PA  (AMBIENT  PRESSURE)  -  KGM  /  MTR-SEC  SOR’D 
ALPHA  (ABSORPTION)  -  1  /  CM 

100  MAMELIST/FIPST/GMA,TA,PAf  ALPHA 
IF  (IFLAG.EQ.l)  ALPHAS. 0 
READ(5, FIRST) 

WRITE  (6,900) 

WRITE(fe, FIRST) 

RH0A=PA/(287.21*TA) 

VFLA=SQRT(GMA*TA*287.21) 

ALL^UNI TS^MUST^BE  ^CONsVsTENT* TO* LNSURE  DIME  NS IONLESS  OUTPUTS • 
THE  UNITS  USED  FOR  THESE  PARTICULAR  RUNS  WERE  AS  FOLLOWS: 

ZMA  (MACH  NUMBER)  -  DI  MENS  IONLFSS  Tf,K1  i  ccc 

JM  (NUMBER  OF  ROWS  IN  THF  Y-DIRECTION)  -  DIMENSIONLESS 

ZL  (TOTAL  LENGTH  OF  INTEREST  IN  X-DTR  ECTI ON)  -  CM 

ZH  (TOTAL  LENGTH  OF  INTEREST  IN  Y-DI RECTIQN )  -  CM 

TAU  (KINETICS)  -  SECONDS 

Z 10  (INTENSITY)  -  WATTS  /  MTR  SOR’D 

SGM  (STANDARD  DEVIATION  OF  GAUSSIAN  BEAM)  -  CM 


200 


READ  (5 ,910, END =999)  ZMA, JM , ZL , ZH  ,TAU , Z I  0, SGM 
IF  (IFLAG.EQ.l)  ZI0=1.0 
YDEL=ZH/( JM-1 ) 

IE(ZMA.LT.l.O)  XDEL=YDEL 
IF  (ZMA.LT.l .0)  GO  TO  200 
ZMU=ARSTN( 1.0/ ZMA) 

XDEL-YDEL/TAN(  ZMU) 

IFlZMA^GTato)  BTA=SQRT(ZMA**2-I  .) 

IF ( ZMA. LT .1.0)  BTA=SQRT( l.-ZMA**2) 

U=ZMA*VEL A  _  „ 

ACH=(GMA-1.0)*YDEL* ALPHA/ 2.0 
ARHO= ZMA/ ( VELA**3*BTA*RH0A) 

AWAK= ( GMA-1. )/( VELA**2*U*RH0A) *XDEL*ALPHA 


onoo  nonon  o  ooo 
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IF  (TAU.LT.i.OE-25)  GO  TO  210 
XPON=XDEL/(U*TAU*10Q) 

IF  (XPON. GT. 35.0)  GO  TO  210 
XL AG= l / EXP (XP 3N) 

GO  TO  220 
210  XLAG=0. 0 
220  C0NTINUF 
I MX= I M 

IM  WILL  BE  INCREMENTED  BY  ONE  TO  INITIALIZE  CALCULATIONS  IN  DIDER 
I M=IM+1 

IF(  JM.GT.JN.OR.IM.GT.IM  GO  TO  240 
GO  TO  250 
240  WRITE  (6,920) 

250  CONTINUE 

CALL  DIDER  ( IN, JN, IFL AG) 

CALL  WRKR  (IN,IMX, JFLAG,KFLAG,LFLAG) 

900  FORMAT! *  1  *  ,T50 , • INPUT  PARAMETERS*,///) 

920  FORMAT( *1*1* All5wED°5imInSIOn|Ff6r  THE  ARRAY  LESS  THAN  REQUIRED* ) 
GO  TO  100 
999  STOP 
END 


SUBROUTINE  DIDER 


SUBROUTINE  DIDER  ( IN,  JN,  I FL  AG  )  _ 

COMMON/ONE/ 1 M,  JM, BTA,ZMA,XDEL,YDEL,ZIO,SGM 
COMMON/TWO/ACH, AWAK,XLAG, ZH, ZL 

DIMENSION  ZI ( IN, JN) ,  IN  AND  JN  CONSISTENT  MATRIX  DIMENSIONS 

AND  Xll(I),  X 12  (  I )  ,  ETA1U),  AND  ETA2CI),  I  EQUAL  TO  THE  NUMBER 
CONSTANT  HEAT  RELEASE  DISTRIBUTIONS  DESIRED. 

DIMENSION  ZI (65, 65  I, XII (8 ) , X  I? ( 8 ) , ET A1 <  8 ) , ETA2( 8 ) , 1  IRC  8 ) 

DIMENSION  CHR(IN,JN),  CHRX(IN,JN),  CHRY(IN,JN),  WAK(INtJN),  ZIC(IN,JN) 
IN  AND  JN  THE  SAME  AS  ZI(IN,JN) 

COMMON/ FOUR/ CHI (65,65) ,CHR X( 65, 65) ,CHRY( 65,65) ,WAK( 65,65) ,ZIC (65,6 


O* 


onnnooo 


lB**?n)*(GMA-1.0V*C/(6.283t8*BTA) 

IF  (IFLAG.EQ.l)  GO  TO  15 

FILL  IN  ENERGY  DENSITY  (GAUSSIAN  HEAT  RELEASE  DISTRIBUTION) 

DO  5  1=1 t IM 
DO  5  J=1,JM 
ZI(I,J>=0.0 
5  CONTINUE 
X  Z IO=l-XLAG 
Xt=ZL/2.0 
Y l=ZH/2  .0 
00  10  I =2  *  I M 

DO  10  J=1,JM 
YY=( J-l )*YDEL 

y y_/ T  _  p 1*X0FL 

s= T  (  x X- X  l )  **  2  f  ( Y  Y- Y 1 )  * * 2 )  /  (  2  •  0*  S  GM* *  2 ) 

IFCS.GT.35.0)  GO  TO  10 
ZI(It J)=ZIO*XZIO/EXP(S) 

10  CONTINUE 
GO  TO  45 

FILL  IN  ENEPGY  DENSITY  (CONSTANT  HFAT  RELEASE  DISTRIBUTION) 

INPUT  ZIR(I)  -  tfATTS  /  MTR  SQR'D-CM  FOR  EACH  CONSTANT  HEAT  RELEASE 
DISTRIBUTION!  ( I) .  ALL  UNITS  MUST  BE  CONSISTENT  TO  INSURE  DIMEN¬ 
SIONLESS  OUTPUTS. 

15  READ  ( 5,400)  IZIR(I),I=1,8) 

DO  20  1=3,23 
DO  20  J=2,17 
ZKI,  J)=ZIR(1) 

20  CONTINUE 

DO  22  1=11,15 
DO  22  J  =9,10 
ZKI,  J)=ZIR(2) 

22  CONTINUE 

DO  24  1=26,46 
DO  24  J  =  2 *  1 7 
ZKI,  J)=ZIR(3) 

24  CONTINUE 

DO  26  1=34,38 
DO  26  J=9 , 10 
ZI(I, J)=ZIR(4) 

26  CONTINUE 

DO  28  1=3,23 


ooo  o  on  oooo 


7 


t 


28 


30 


32 


34 


45 


50 


00  28  J=20,35 
ZKIt  J)=ZIR<5) 

CONTINUE 
DO  30  1=11,15 
DO  30  J=27,28 
ZI  (I,J)=ZIR(6) 

CONTINUE 
DO  32  1=26,46 
DO  32  J  =  20, 35 
ZKIt  J1'=ZIR(7) 

CONTINUE 
DO  34  1=34,38 
DO  34  J=27,28 
ZI ( I , J)  =  Z I k{ 8 1 
CONTINUE 

ZEROIZE  COMPUTATIONAL  MATRICES 
KINETICS. 

DO  50  J=1 , JM 
ZlCUt  J1=0.0 
WAKC1, J)=0.0 
00  50  I=2,IM 

zicu,  j)=o.o 

WAK(  ,  J)  =  0.0 
CHR(I, J)=0.0 
CHRX( I , J) =0. 0 

CHRY  ( I ,  J )  =0.  0  ^  it 

ZIC(I, J)=Zir(I-l, J)*XLAG  +  ZI(I,J1  _  ,  i%i/o 
WAK(I  ,J)=WAK(I-l,  J)  +  CZIC(I,JKZIC(I-1,  JK/2 
CONTINUE 

COMPUTE  WAK*=  TERM  OF  DENSITY  PERTURBATION. 


AND  CALCULATE  HEAT  RELEASE  FROM 


DO  60  J=1 , JM 
DO  60  1=2, IM 

60  WAK( I ,J)=AWAK*WAK( I ,  J  ) 
IF(IFLAG.EO.l)  GO  TO  200 
IFUMA.LT. 1.0)  GO  TO  110 


CALCULATE  CONTRIBUTIONS  FROM  CHARACTERISTICS  FOR  SUPERSONIC  FLOW 


DO  100  I=2,IM 
00  100  J=1 , JM 
JL  =  J 
JF=J 
L  =  I 

DO  80  N=2 , 1 


noooonnoonoonoo 


c  ft 


L  =  L-1 
JR= JR»l 

CHRX(  !  »  J)  =CHRK  (  I?  J  )♦  ACH*(  Z I C ( L  ,  JR  )*Z  IC  ( L  ♦  1 ,  JR  +  1 ))  /2  •  0 
70  JL= JL+1  m  _ 

CHRY(  It  J)=CHRYf  1  ,  J)+ACH*(ZIC(L,  JL)+ZIC<  L  +  l ,  JL-1  ) )  /2.0 
80  CONTINUE 

CHR( I  *  J)=CHRX( I  *  JUCHRYI It J) 

CHPYU  ,  J)  =CHRY(  f  ,  J)-CHRX(  I ,  J) 

CHRXII,  J)=-CHR(I,J) 

100  CONTINUE 
RETURN 

CALCULATE  CONTRIBUTIONS  FROM  SOURCE  TERMS  FOR  SUBSONIC  FLOW, 

110  00  150  1=2, IM 
DO  150  J=1 , JM 
I I 1=MAX0( 2,1-15) 

II2=MIN0(IM,U15> 

DO  150  11=111,112 

DO  150  J J=1 , JM  %%  _  ^ 

I F < ( 1 1 • EQ • I > . AND. ( J J • FQ. J ) )  GO  TO  150 
UX=(I-II)*XDEL 

CHfiWt  J)=CHR(  I?J)-C1(UX,UY)*ZIC< II,  J.J1 

CHRYUl  jl^CHRY?  1 1  jltCl (UY,UX)*ZIC( II  ,  J  J 1 
150  CONTINUE 
RETURN 

200  CONTINUF 


CALCULATE  CONTRIBUTIONS  FROM  SOURCE  AMD  WAKE  TERMS  FOR  SUBSONIC 
FLOW  WITH  CONSTANT  HEAT  RELEASE  DISTRIBUTION. 

INPUT  THE  ORDINATES  OF  EACH  CONSTANT  HEAT  RFLEASE  DISTRIBUTION: 
XIKI)  AND  X 12  (  I )  IN  THE  X-DIRECTION  - ;  CM. 

ETAim  AND  ETA2II)  IN  THE  Y-DIRFCTION  ~  CM. 


NOTE*  A. 


B. 


IN 

1. 

2. 
IN 
1. 
2. 
IN 
1. 


THF  X  (FLOW)  DIRECTION: 

MINIMUM  VALUE  OF  XIKI)  IS  0.0. 

MAXIMUM  VALUE  OF  XI2CI)  IS  ZL.  ^ilTpni  THC. 

THE  Y-DIRECTION  (HEAT  DISTRIBUTION  ABOVE  CENTERLINF) 
MI  MI  MUM  VALUE  OF  ETAl(I)  IS  0.0. 

MAXIMUM  VALUE  OF  ETA2U)  IS  ZH/2.  Brinll  „llTeill 
THE  Y-DI RECTION  (HEAT  DISTRIBUTION  BELOW  CENTERLINE) 
MTNTMUM  VALUE  OF  ETAl(I)  IS  -ZH/2. 


ononn 


k 


IN  THF^Y-OI REcVVoN^fHEAT'DISTR  I  BUT I  ON  ACROSS  CENTERLINE ) 

i:  km  »E  of  mn\  if M2.* 


ETA2U 


ZH/2. 


READ  (5,410)  <XIl(I),XI2m,ETAl<ntFTA2<n,I=l,8) 

DO  300  K=l,8 
DO  250  1=2, IM 
DO  250  J=l » JM 
UX=XDEL*( 1-2 ) 

UY=(ZH/2.0-YDEL*( J-l  ) ) 

DO  250  I MF=1 , 5 

r  ljd  |  r  ^  i  sfHR  (I«J1+F((  UX-X  I2(K))»(  UY-  (  - 1 )  **  I  MG*FT  A 1  ( K  )  * 1  MG*ZH )  *BT  A  , 
1  Z IR( K ) )  ~F(  ( UX-X I  * ( K) ) ? ( UY- ( -l ) ** IMG*ET  A 2  ( K )  +1  MG*  ZH)*BT A , ZIR ( K) |-F 
2 ( (UX-XI 1 (K) ) t  (  UY- f -1 ) i*IMG*ET A1 ( K )  +  I MG*Z H )*BT  A, Z IR ( K ) ) +F((UX-XI1(K 

1 ) «ZIR(K))-F( (JY-(-l) **l MG*  FT A 1 ( K ) +1 MG*ZH )*BTA*(UX  — X 1 2  ( K  ) ),ZIR(K))— 
2F I (UY-(-l) **IMG*5TA2( K)tIMG*ZH)#BT  A,(UX-XT1 ( K I > *  Z ! R<  K) )  +F( (UY- (-1 ) 
3**IMG*ETA1( K I HM3*ZH)*BTA,(UX-XI1(K)),ZIR(K)) 

250  CONTINUE 
300  CONTINUE 

DO  350  1=2, IM 
DO  350  J=1 , JM 
CHRX(I,J)  =-CHR  ( I  ,  J  ) 

^50  CONTINUE 

400  FORMAT  (4E12.2) 

410  FORMAT  (4F10.6) 

RETURN 

END 


SUBROUTINE  WRK* 


SUBROUTINE  WRKR  ( IN, I  MX, JFL AG, KFL AG,LFL AG) 

COMMON/ONE/ IM, JM , BT A, ZMA, XDEL, YDEL, ZIO, SGM 
COMMON/THREE /GMA , TA, P A,U, ARHO, T AU 

DIMENSION  CHR  ( IN , JN) «  CHRX(IN,JN),  CHR Y ( IN , JN) ,  WAK (  IN , JN )  ,  ZIC(IN,JN) 
IN  AND  JN  THE  SAME  AS  ZI(TN,JN)  IN  THE  MAIN  PROGRAM 

COMMON/ FOUR/CHRC  6  5  ,65  )  , CHRX C  65 , 6 5  )  ,CHP  Y  (  65 , 65 )  , W  AK (  6 5,  6 5  )  ,  Z1C(65,6 
15) 


ooo  non 


*  ■  x 


SEE  LISTING  OF  SUBROUTINE  CONTUR  FOR  INFORMATION  ON  TITLE(12),  LTG( 3 ) * 
CL( 8 )  AND  OTHER  PERTINENT  PARAMETERS. 


IF  (IFLAG.FQ.l)  GO  TO  10 


RFAL*8 
10NT0U 
20NTOU 
3CITY 
4CITY 
5HIXX 
6YY  • 
7SSI A* 
GO  TO 
10  REAL* 8 
1 ONTOU 
20NT0U 
3  CITY 
4CITY 
5HIXX 
6YY  • 

7 ST  AN" 


MACH 


•/ 

«/ 


TITLEt 12) 2TLINE1(6,6)/ 
,  1 RS  S 

,'RS  *, 

, 'CONTOURS*  , 
t  *  CONTOURS  *  , 

\\ 

•N  '  ,* 

20 

TITLE (12). 

,  *  RS  * 

t  *RS  * 

,  *  CONTOURS* 

, 'CONTOURS* 

t  *  * 

•  )  *  , 

•T  '  - 


TL INE2 ( 5  )  / 


I* 

IS* 

PERTUR* 

PERTUR* 


, ' SODENSIT 
t ' OPR  ESSUR 
,  *  BAT I ON  V 
,*BATION  V 


TL INE1 (6,6)/ 

, 

, 

, 

, 


TLINE2I5)/ 


50 


100 


'/ 

MACH  * / 

20  L0GICAL*1  LTG<3)/. TRUE., .TRUE. ..TRUE./ 

REAL*4  CL ( 16) /-60. 0,-55. 0,-50. 6, -45. 0,-40. 

1 0, -1 5, 0,-10. 0,-5. 0,0.0, 5.0 t 10.0,15.0/ 

JL= JM/2 
JLi  =  (  JM4-1 )  /2 

APRES=l .0E+06*ZMA**2/BTA/U/PA 
AYVFL=1 .0E+06/GMA/U/PA 
AXVEL= AYVEL /BT  A 
DL=YDEL /XDEL 
CP=GMA*159.56/(GMA-1) 

TO=TA*(  l+(GMA-l)*ZMA**2/2.  ) 

I MXX= I M-2 
JMX= J  M- 1 
DO  50  1*1,5 
J  =  U6 

TITLE( J)=TLINE2( I  ) 

CALL  REREAD 
READ  ( 5.950)  AMA 
RFAD  (99,960)  TITLE(12) 

ZTC  GIVES  DENSITY  PERTURBATIONS 

DO  100  1=2, IM 
DO  100  J=l , JM 

ZIC(I-1,J)=1.DE+06*(ARH0*CHR(I,J)-WAK(T,J) ) 
IF  (IFLAG.EQ.l)  JL1  =  JM 
IE  ( JFLAG.EQ.O)  GO  TO  160 


*  . 

•  ,*  ' , 

HEAT*,*  INPUT  -• , 

*  PHI  Y 

«  GAU 

I • , • SODENSIT* 

, « Y  C 

IS*,*  OPRESSUR ' 

,*E  C 

PERTUR* , • RATION  V* 

,'ELO 

PERTUR* , ' BATION  V* 

,'ELO 

*,'  * 

,  *  P 

•  *  •  , 

HEAT*,*  INPUT 

'  PHI  Y 
•  CON 

,-35.0,-30.0,-25.0 

,-20. 

,'Y  C 
, » E  C 

t*FLO 
t*  FLO 


oonno 


to 


110 

115 

120 


130 

140 

150 


IF  (JFLAG.EQ.2)  GO  TO  140 
IF(ZMA.LT.l.O)  SO  TO  110 
WRITE  ( 6*900 )  ZMA 
GO  TO  120 

IF  (IFLAG.EQ.l)  GO  TO  115 
WRITE  (6,901)  ZMA 
120 

(6,902) 

(6,910) 

(6,920) 

( 6,930) 

1*1, IMX 


ZMA 
TAU  . 

I  MX , JM 


GO  TO 
WRITE 
WRITE 
WRITE 
WRITE 
00  130 

WRITE  (6,940)  (ZICU, 

IF  (JFLAG.EQ.l)  GO  TO 
DO  150  1*1,6 
TITLF( I )=TLINEl ( I ,1) 

CALL  CONTUR  ( Z I C , I  MX , JM, I N, C L, -8 , TITLE, 6 , 


J)  , J=1 , JL1) 
160 


160 


200 


210 

215 

220 


222 

224 

230 


235 


ZIC  GIVES  PRESSURE  PERTURBATIONS 
CHR  GIVES  X-VELOCITY  PERTURBATION 
WAX  GIVES  Y-VELOCITY  PERTURBATION 


DO  200  1=2, IM 
DO  200  J=1 , JM 
ZIC(I-1,J)=APRES*CHR(I ,  J) 

C  HR ( I - 1 , J)=CHRX( I, J)*AXVEL 
WAKd-l  ,  J)=CHRY( I, J)*AYVEL 
CONTINUE 

IF  (KFLAG.EQ.3)  GO  TO  260 
IF  (KFLAG.E0.2)  GO  TG  245 
IF  ( JFLAG.EQ.l. GR.JFLAG.EQ. 3) 
IF  (ZMA.LT.l.O)  GO  TO  210 
WRITE  (6,900)  ZMA 
GO  TO  220 

IF  (IFLAG.EQ.l)  GO  TO  215 
WRITE  ( 6,901 )  ZMA 
GO  TO  220 

ZMA 
TAU 

IMX, JM 


GO  TO  222 


WRITE 
WRITE 
WRITE 
WRITE 
GO  TO 
WRITE 
DO  230 


(6,90?) 
(6,910) 
(6,920) 
(6,931  ) 
224 

(  6,932) 

1  =  1, IMX 


WRITE  (6,940)  ( Z IC ( I , J ) , J=1 , JL1 ) 
WRITE  (6,933) 

DO  235  1=1, IMX 

WRITE  (6,940)  (  CHR( I , J) , J  =  1 , JL1 ) 


6  ,  LT  G ) 


2  40 

245 

250 


252 

254 

260 


300 


400 


600 

605 

610 


612 

614 

620 


DO  240  1*1 f IMX 

WRITE  (6,940)  ( W AK( I , J) , J  =  1 , JL1 ) 

IF  (KFLAG.EQ.l)  GO  TO  260 
DO  250  1  =  1,6 

TITLE(I)  =  TLTNE1(  1,2)  „  ,  #  _ 

CALL  CONTUR  ( Z IC , I  MX, JM, I N, CL, -8 , T ITLE , 6 ,6 , LTG) 

DO  252  1=1,6 

TITLE ( I )=TLINE 1 ( I , 3) 

CALL  CONTUR  ( CHR , IMX , JM, I N , C L , -8 , TITLE , 6 ,6 , LTG) 

DO  254  1*1,6 
TITLE( 11  =  TL I NF1(  1,4) 

CALL  CONTUR  ( WAK , f MX , JM, I N, C L, -8 , TITLE , 6 , 6 , LTG) 

IF  (LFLAG.EQ.O)  RETURN 
DO  300  J=1 , JM 

CHRX( 1 , J) =  C  HP  (  2  ,  J  ) -CHR ( 1 , J ) 

CHRX(IMX»J)=  IMX,  J)-CH&  (IMXX,  J) 

DO  300  I =2 , I MXX 

CHPX( I, J)=(CHR(  I  *1, J)-CHR( 1-1, J) )/2.0 

CONTINUE 

DO  400  1  =  1  §  I  MX 

CHRY(I,1)=(WAK( I ,1  )-WAK( I, 2))/DL 
CHRY(  I,  JM  )=(  WAK(  I  ,  JMX  ) -WAK  ( I  ,  JM)  ) /DL 
DO  400  J=2  *  JMX 

CHRY ( I ,  J )  =  f W AK( I , J-l )  -WAK ( I , Jf 1 )  )/(2.0*DL) 

CONTINUE 

CHRX  GIVES  PHI  XX 
CHRY  GIVES  PHI YY 

IF(JFLAG?EQ?i!oR?JFLAG,EQ.3.0R.KFLAG,EQ.l,OR.KFLAG«EQ.3)  GO  TO  612 
IF  (ZMA.LT.1.0)  GO  TO  600 
WRITE  (6,900)  ZMA 
GO  TO  610 
IF  (IFLAG.FQ.l) 


GO  TO  605 


WRITE 
GO  TO 
WRITE 
WRITE 
WRITE 
WRITE 
GO  TO 
WRITE 
DO  620 


(  6,901  ) 
610 

(  6,902) 
(6,910) 

(  6,920) 

(6.935) 
614 

( 6.936) 
1=1, IMX 


ZMA 

ZMA 

TAU 

IMX, 


JM 


WRITE  (6,940)  ( C HRX( I , J ) , J  =  1 , JL 1  ) 
WRITE  (6,937) 

DO  630  1=1, IMX 
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630  WRITE  (6,940)  (CHRY(I »J) » J=1 , JL1 ) 

IF  (LFLAG.EQ.l)  RFTURN 
635  00  640  1=1,6 

640  TITLE(  I)=TUNE1(  I  ,5)  _  _ _ _ .  _  ,  .  ,  Trl 

CALL  CONTUR  ( CHRX , IMX , JM, I N ,CL , -8 , TITLE , 6, 6 ,LTG) 

DO  650  1=1,6 

650  TITLE(I)  =  TLIN51U  ,6)  TTT,  _  ,  .  .  Tr% 

900  GAUSSIAN  HEAT 

9011FORMATE(?lI>!//lox!;SUBSnNIC  CASE  (M=*,F7.5,*I  AND  GAUSSIAN  HEAT  RE 

DOZ^ORMa/oIn/ZOOX^SUBSONIC  CASE  <M=*,F7.5,M  AND  CONSTANT  HEAT  RE 

9  1 0^  FORMA  T°  (  30x1  •'xiNF  Til?  S^ !  TAU=  •»F10.3»*J  AND  NO  WALLS*.  US 
-  ( 30X, •MATRIX:  *,I3,*  COLUMNS  AND  *, I */ /I 


9 20  FORMAT 

930  FORMAT 

931  FORMAT 
93?  FORMAT 

1//) 

933  FORMAT 
Iff ) 

934  FORMAT 
Iff) 

935  FORMAT 

936  FORMAT 

937  FORMAT 
940  FORMAT 
950  FORMAT 
960  FORMAT 

RFTURN 

END 


(  30X, * KINETIC  S  ( TAU= 1 ,E 10 . 3 , • )  AND  NO  WALLS*,//) 

(30x5* MATRIX:  *,I3,*  COLUMNS  AND  *,I3T*  R°WS*,//) 

(50x1  'TOTAL  DENSITY  CHANGFS,  TOP  HALF  OF  LASER  BEAM*,//) 

<  46X,  *  PRESSURE  _  P.FRTURBATI  ONS,  TOP  H£kF  Of  LASER  BEAM  •  , // 


( « 1 1 ,45X, • PRFSSURE  PERTURBATIONS,  TOP  HALF 

( *  1*  *45X, *U  PERTURBATION  VELOCITY,  TOP  HALF 

(*1*,45X,*V  PERTURBATION  VELOCITY,  TOP  HALF 

(66X,  *  PHI XX  *  »//  )  , 

(*l*,65X, 1  PHI  XX  * , //) 

(*1*,65X,* PHI YY* , // ) 

(1X,1P13E10.2) 

( F8  *6 ) 

(A8) 


OF  LASER  BEAM* 
OF  LASER  BEAM* 
OF  LASER  BEAM* 
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